<jupyter_start><jupyter_text>Stifness matrix of the P2 triangle element.<jupyter_code>x,y=var("x,y")<jupyter_output><empty_output><jupyter_text>The barycentric coordinates on the reference triangle:<jupyter_code>lamb=[1-x-y, x,y]<jupyter_output><empty_output><jupyter_text>The P2 basis functions:<jupyter_code>phi=[2*l*(l-1/2) for l in lamb]+[4*lamb[i]*lamb[(i+1)%3] for i in range(0,3)]<jupyter_output><empty_output><jupyter_text>Verify that basis functions value is zero or one where it should be:<jupyter_code>points=[[0,0],[1,0],[0,1],[1/2,0],[1/2,1/2],[0,1/2]]
A=matrix(6,6)
for i in range(0,6):
    f=phi[i]
    for j in range(0,6):
        A[i,j]=f(x=points[j][0],y=points[j][1])
assert(A==identity_matrix(6,6))<jupyter_output><empty_output><jupyter_text>Gradients of  basis functions on the reference triangle:<jupyter_code>grads=[[diff(f,x),diff(f,y)] for f in phi]<jupyter_output><empty_output><jupyter_text>Gradients at quadrature points (midpoints of the edges), on the reference triangle:$gq_{i,j}$= gradient of ith function at jth quadrature point.<jupyter_code>quad=[[1/2,0],[1/2,1/2],[0,1/2]]
gq=[[(f[0](x=q[0],y=q[1]),f[1](x=q[0],y=q[1])) for q in quad]  for f in grads]<jupyter_output><empty_output><jupyter_text>Transformation matrix J from reference to current triangle, and its inverse:<jupyter_code>x0,x1,x2,y0,y1,y2=var("x0,x1,x2,y0,y1,y2")# current triangle
J=matrix([[x1-x0,x2-x0],[y1-y0,y2-y0]])
det=J.determinant().simplify_rational()
Jinv=J.inverse().simplify_rational()
JinvDet=Jinv*det #factorize determinant.
# transpose of J^-1 multiplied by determinant:
JinvDetTrans=JinvDet.transpose()<jupyter_output><empty_output><jupyter_text>Verify that 'everything' is ok on the reference triangle:<jupyter_code>assert(det(x0=0,x1=1,y0=0,y1=0,x2=0,y2=1)==1)
assert(JinvDetTrans(x0=0,x1=1,y0=0,y1=0,x2=0,y2=1)==identity_matrix(2,2))<jupyter_output><empty_output><jupyter_text>Exact integration of a function f(x,y) on the reference triangle:<jupyter_code>integreftr=lambda f: f.integral(x,0,1-y).integral(y,0,1)<jupyter_output><empty_output><jupyter_text>Stiffness matrix on reference triangle, exactly computed:<jupyter_code>M=matrix(QQ,6,6)
for i in range(0,6):
    for j in range(0,6):
        gigj= grads[i][0]*grads[j][0]+grads[i][1]*grads[j][1]
        M[i,j]=integreftr(gigj)<jupyter_output><empty_output><jupyter_text>Stiffness matrix on any triangle, numerically computed:<jupyter_code>MN=matrix(SR,6,6)
Ji=JinvDetTrans
for i in range(0,6):
    for j in range(0,6):
        s=0
        for k in range(0,3):
            V=Ji*vector(gq[i][k])
            W=Ji*vector(gq[j][k])
            s+=V.dot_product(W)
        MN[i,j]=s
MN*=(1/(6*det))<jupyter_output><empty_output><jupyter_text>Verify that stiffness matrix as computed by quadrature is equal to the  'exact' matrix on the reference element. Recall that our numerical quadrature is exact for 2nd order polynomial and that, here, we compute in $\mathbb{Q}$ (all computations are exact).<jupyter_code>assert(MN(x0=0,x1=1,y0=0,y1=0,x2=0,y2=1)==M)<jupyter_output><empty_output><jupyter_text>Verify also that if v=[1,1,1,1,1,1], MN.V==0 (on any triangle):<jupyter_code>x=vector([1,1,1,1,1,1])
assert((MN*x).simplify_full().is_zero())<jupyter_output><empty_output><jupyter_text>How is the stiffness matrix on the reference triangle?<jupyter_code>print(M)<jupyter_output>[   1  1/6  1/6 -2/3    0 -2/3]
[ 1/6  1/2    0 -2/3    0    0]
[ 1/6    0  1/2    0    0 -2/3]
[-2/3 -2/3    0  8/3 -4/3    0]
[   0    0    0 -4/3  8/3 -4/3]
[-2/3    0 -2/3    0 -4/3  8/3]<jupyter_text>Define a function that will return the stiffness matrix, for x=[,,], y=[,,] containing the coordinates  of the summits of a triangle:<jupyter_code>def Stiff(x,y):
    return MN(x0=x[0],x1=x[1],y0=y[0],y1=y[1],x2=x[2],y2=y[2])<jupyter_output><empty_output><jupyter_text>Do we again compute the same matrix on the reference triangle?<jupyter_code>x=vector([0,1,0])
y=vector([0,0,1])
assert(Stiff(x,y)==M)<jupyter_output><empty_output><jupyter_text>The Stiffness matrix should be invariant under dilatations and rotations. Verify this, for some values, dilating and rotating the reference element.We must rotate by angles whose sinus and cosinus are Algebraic Numbers, and compute in (real) Algebraic Numbers (AA for SageMath), so as to make again exact computations (and we choose dilatation factor in $\mathbb{Q}$ or in Algebraic numbers):<jupyter_code>tetas=[pi/12,pi/6,pi/4,pi/3,pi/2]
for teta in tetas:
    Rot=matrix(AA,[[cos(teta),-sin(teta)],[sin(teta),cos(teta)]])
    for dilat in [1/10,3/10,7/10,11/10,sqrt(2),2,sqrt(5),4,6,71/10]:
        RotD=dilat*Rot
        V=[RotD*vector(AA,[x[i],y[i]]) for i in range(0,3)]
        assert(Stiff([s[0] for s in V],[s[1] for s in V])==M)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.17156040268456377, 'fraction_numerical': 0.05432046979865772, 'mean_word_length': 6.214826021180031, 'pattern_counts': {'":': 0, '<': 68, '<?xml version=': 0, '>': 68, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '17', 'dir': 'jupyter-structured-clean-dedup', 'id': '402996', 'license': 'MIT', 'n_tokens_mistral': 2085, 'n_tokens_neox': 2016, 'n_words': 384, 'path': 'FeStiff/sage/Stiff.ipynb', 'repo_name': 'wolfv/BenchmarksPythonJuliaAndCo', 'size': '4768'}