Tähän liittyviä skannattuja kalvoja

Kolmio1
Kolmio2
Kolmio3
Kolmio4

>

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]]);

cartesian := vector([1, x, y])

natural := vector([L1, L2, L3])

transform := matrix([[1, 1, 1], [x1, x2, x3], [y1, ...

> op(cartesian)=evalm(transform&*natural);

vector([1, x, y]) = vector([L1+L2+L3, x1*L1+x2*L2+x...

> matrix(3,1,lhs(%))=matrix(3,1,rhs(%));

matrix([[1], [x], [y]]) = matrix([[L1+L2+L3], [x1*L...

> itr:=inverse(transform):

> LL:=evalm(itr&*cartesian);

LL := vector([-(x2*y3-x3*y2)/(-x2*y3+x3*y2+x1*y3-x1...
LL := vector([-(x2*y3-x3*y2)/(-x2*y3+x3*y2+x1*y3-x1...
LL := vector([-(x2*y3-x3*y2)/(-x2*y3+x3*y2+x1*y3-x1...
LL := vector([-(x2*y3-x3*y2)/(-x2*y3+x3*y2+x1*y3-x1...
LL := vector([-(x2*y3-x3*y2)/(-x2*y3+x3*y2+x1*y3-x1...

> twoA1:=det(matrix(3,3,[1,1,1,x,x2,x3,y,y2,y3]));

twoA1 := x2*y3-x3*y2-x*y3+x*y2+y*x3-y*x2

> twoA:=det(matrix(3,3,[1,1,1,x1,x2,x3,y1,y2,y3]));

twoA := x2*y3-x3*y2-x1*y3+x1*y2+y1*x3-y1*x2

> simplify(twoA1/twoA-LL[1]);

0

This was a proof that the "natural coordinates" and the "area coordinates" agree.

> A:='A':LL[1];

(-x2*y3+x3*y2)/(-x2*y3+x3*y2+x1*y3-x1*y2-y1*x3+y1*x...
(-x2*y3+x3*y2)/(-x2*y3+x3*y2+x1*y3-x1*y2-y1*x3+y1*x...

> 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]);

L1 := 1/2*(-x2*y3+x3*y2)/A+1/2*(y3-y2)*x/A-1/2*(x3-...

L2 := 1/2*(x1*y3-y1*x3)/A-1/2*(y3-y1)*x/A-1/2*(-x3+...

L3 := -1/2*(x1*y2-y1*x2)/A+1/2*(y2-y1)*x/A+1/2*(-x2...

> abc1:=coeffs(L1,[x,y]);abc2:=coeffs(L2,[x,y]);abc3:=coeffs(L3,[x,y]);

abc1 := 1/2*(-x2*y3+x3*y2)/A, -1/2*(x3-x2)/A, 1/2*(...

abc2 := 1/2*(x1*y3-y1*x3)/A, -1/2*(-x3+x1)/A, -1/2*...

abc3 := -1/2*(x1*y2-y1*x2)/A, 1/2*(-x2+x1)/A, 1/2*(...

> 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);

1/2*(-x2*y3+x3*y2)/A+1/2*(y3-y2)*x1/A-1/2*(x3-x2)*y...

> simplify(subs(x=x1,y=y1,LL[1]));subs(x=x2,y=y2,twoA1/twoA);subs(x=x3,y=y3,twoA1/twoA);

1

0

0

Hence L1(x1,y1)=1 and is 0 on other corners, similarly for L2 and L3, ie. L[i](a^j) = delta[i,j]

(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];

1-x-y

> 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;

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])));

matrix([[b1^2+c1^2, b1*b2+c1*c2, b1*b3+c1*c3], [b1*...

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])));

M := A*matrix([[b1^2+c1^2, b1*b2+c1*c2, b1*b3+c1*c3...

We have all these quantities available in terms of the coordinates of the vertices.

> A=twoA/2;

A = 1/2*x2*y3-1/2*x3*y2-1/2*x1*y3+1/2*x1*y2+1/2*y1*...

> abc1;abc2;abc3;

1/2*(-x2*y3+x3*y2)/A, -1/2*(x3-x2)/A, 1/2*(y3-y2)/A...

1/2*(x1*y3-y1*x3)/A, -1/2*(-x3+x1)/A, -1/2*(y3-y1)/...

-1/2*(x1*y2-y1*x2)/A, 1/2*(-x2+x1)/A, 1/2*(y2-y1)/A...

> L[1];

a1+b1*x+c1*y

> L1:='L1';L2:='L2':L3:='L3';

L1 := 'L1'

L3 := 'L3'

> op(transform);

matrix([[1, 1, 1], [x1, x2, x3], [y1, y2, y3]])

> onexy:=evalm(transform&*[L1,L2,L3]);x:=simplify(subs(L3=1-L2-L1,onexy[2]));y:=subs(L3=1-L2-L1,onexy[3]);

onexy := vector([L1+L2+L3, x1*L1+x2*L2+x3*L3, y1*L1...

x := x1*L1+x2*L2+x3-x3*L2-x3*L1

y := y1*L1+y2*L2+y3*(1-L2-L1)

> J:=jacobian([x,y],[L1,L2]);det(J);twoA;

J := matrix([[-x3+x1, x2-x3], [y1-y3, y2-y3]])

x2*y3-x3*y2-x1*y3+x1*y2+y1*x3-y1*x2

x2*y3-x3*y2-x1*y3+x1*y2+y1*x3-y1*x2

Hence dx*dy = 2*A*dL1*dL2

> 2*'A'*int(int(Li*Lj,Lj=0..1-Li),Li=0..1);

1/12*A

> 2*'A'*int(int(Li*Li,Lj=0..1-Li),Li=0..1);

1/6*A

> M:=(A/6)*matrix(3,3,[1,1/2,1/2,1/2,1,1/2,1/2,1/2,1]);

M := 1/6*A*matrix([[1, 1/2, 1/2], [1/2, 1, 1/2], [1...

This is the local "mass matrix".

Quadratic shape functions

> LL[1]:

> L[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[0]

psi[2] := 4*L[0]*L[1]

psi[3] := 4*L[1]*L[2]

> psi[1]:=4*L[2]*L[3];psi[2]:=4*L[1]*L[2];psi[3]:=4*L[1]*L[2];

psi[1] := 4*L[2]*L[3]

psi[2] := 4*L[1]*L[2]

psi[3] := 4*L[1]*L[2]

>