>
Warning, the name changecoords has been redefined
Area coordinates for a triangle
> restart:with(linalg):with(plots):
Warning, new definition for norm
Warning, new definition for trace
> cartesian:=vector([1,x,y]);natural:=vector([L1,L2,L3]);transform:=array([[1,1,1],[x1,x2,x3],[y1,y2,y3]]);
> op(cartesian)=evalm(transform&*natural);
> matrix(3,1,lhs(%))=matrix(3,1,rhs(%));
> itr:=inverse(transform):
> LL:=evalm(itr&*cartesian);
> twoA1:=det(matrix(3,3,[1,1,1,x,x2,x3,y,y2,y3]));
> twoA:=det(matrix(3,3,[1,1,1,x1,x2,x3,y1,y2,y3]));
> simplify(twoA1/twoA-LL[1]);
This was a proof that the "natural coordinates" and the "area coordinates" agree.
> A:='A':LL[1];
>
L1:=subs(denom(LL[1])=2*A,LL[1]);L2:=subs(denom(LL[2])=2*A,LL[2]);
L3:=subs(denom(LL[3])=2*A,LL[3]);
> abc1:=coeffs(L1,[x,y]);abc2:=coeffs(L2,[x,y]);abc3:=coeffs(L3,[x,y]);
> x1:=0:y1:=0:x2:=1:y2:=0:x3:=0:y3:=1:
> plot([[x1,y1],[x2,y2],[x3,y3]]):
> x1:='x1':y1:='y1':x2:='x2':y2:='y2':x3:='x3':y3:='y3':
> subs(x=x1,y=y1,L1);
> simplify(subs(x=x1,y=y1,LL[1]));subs(x=x2,y=y2,twoA1/twoA);subs(x=x3,y=y3,twoA1/twoA);
Hence L1(x1,y1)=1 and is 0 on other corners, similarly for L2 and L3, ie.
(We used just for fun the two different forms of the L-functions.)
> x1:=0:y1:=0:x2:=1:y2:=0:x3:=0:y3:=1:LL[1];
> plot3d(LL[1],x=0..1,y=0..1-x,axes=boxed,orientation=[-120,60]);
> x1:='x1':y1:='y1':x2:='x2':y2:='y2':x3:='x3':y3:='y3':
>
Element stiffness matrix, linear shape functions
> L[1]:=a1+b1*x+c1*y;L[2]:=a2+b2*x+c2*y;L[3]:=a3+b3*x+c3*y;
> matrix(3,3,(i,j)->innerprod(grad(L[i],[x,y]),grad(L[j],[x,y])));
The integration is easy, since the gradients are constants. Thus we just need to multiply this
by the area of our triangle
> A:='A':MS:=A*matrix(3,3,(i,j)->innerprod(grad(L[i],[x,y]),grad(L[j],[x,y])));
We have all these quantities available in terms of the coordinates of the vertices.
> A=twoA/2;
> abc1;abc2;abc3;
> L[1];
> L1:='L1';L2:='L2':L3:='L3';
> op(transform);
> onexy:=evalm(transform&*[L1,L2,L3]);x:=simplify(subs(L3=1-L2-L1,onexy[2]));y:=subs(L3=1-L2-L1,onexy[3]);
> J:=jacobian([x,y],[L1,L2]);det(J);twoA;
Hence
> 2*'A'*int(int(Li*Lj,Lj=0..1-Li),Li=0..1);
> 2*'A'*int(int(Li*Li,Lj=0..1-Li),Li=0..1);
> M:=(A/6)*matrix(3,3,[1,1/2,1/2,1/2,1,1/2,1/2,1/2,1]);
This is the local "mass matrix".
Quadratic shape functions
> LL[1]:
> L[1];
> for i to 3 do phi[i]:=2*L[i]*(L[i]-1/2): od:
> for i to 3 do psi[i]:=4*L[irem(i+1,3)]*L[irem(i+2,3)]:od;
> psi[1]:=4*L[2]*L[3];psi[2]:=4*L[1]*L[2];psi[3]:=4*L[1]*L[2];
>