The problem to solve

by Nasser M. Abbasi  (oct 2009)

In this solution, I start directly by solving for the vector field {u,v} and starting from the general degrees of freedom, and from it by matrix inversion, find the shape function matrix N (in terms of nodal degrees of freedom). This involves inversting a 6 by 6 matrix. But Ok, I am using a computer. By hand, I would use the method I showed in the analytical note part of this assignment which involves inverting only a 3 by 3 matrix.

Needs["Notation`"]
nNodes=3;
nDegreeOfFreedom=6;
Symbolize[(NotationTemplateTag[(Subscript[u, 1])])];Symbolize[(NotationTemplateTag[(Subscript[u, 2])])];Symbolize[(NotationTemplateTag[(Subscript[u, 3])])];Symbolize[(NotationTemplateTag[(Subscript[v, 1])])];Symbolize[(NotationTemplateTag[(Subscript[v, 2])])]
Symbolize[(NotationTemplateTag[(Subscript[v, 3])])];Symbolize[(NotationTemplateTag[(Subscript[a, 1])])];Symbolize[(NotationTemplateTag[(Subscript[a, 2])])];Symbolize[(NotationTemplateTag[(Subscript[a, 3])])];Symbolize[(NotationTemplateTag[(Subscript[a, 4])])]
Symbolize[(NotationTemplateTag[(Subscript[a, 5])])];Symbolize[(NotationTemplateTag[(Subscript[a, 6])])];Symbolize[(NotationTemplateTag[(Subscript[x, 1])])];Symbolize[(NotationTemplateTag[(Subscript[x, 2])])];Symbolize[(NotationTemplateTag[(Subscript[x, 3])])]
Symbolize[(NotationTemplateTag[(Subscript[y, 1])])];Symbolize[(NotationTemplateTag[(Subscript[y, 2])])];Symbolize[(NotationTemplateTag[(Subscript[y, 3])])]

Start by defining the u and v trial functions (linear polynomials in x and y)

nodalDegreesOfFreedom={(Subscript[u, 1]),(Subscript[v, 1]),(Subscript[u, 2]),(Subscript[v, 2]),(Subscript[u, 3]),(Subscript[v, 3])};
generalDegreesOfFreedom={(Subscript[a, 1]),(Subscript[a, 2]),(Subscript[a, 3]),(Subscript[a, 4]),(Subscript[a, 5]),(Subscript[a, 6])};
u=(Subscript[a, 1])+(Subscript[a, 2])x+(Subscript[a, 3])y;
v=(Subscript[a, 4])+(Subscript[a, 5])x+(Subscript[a, 6])y;

set up the  u=X a equation

{b,xMat}=Normal[CoefficientArrays[{u,v},generalDegreesOfFreedom]];
Print[ToString[MatrixForm[{{"u"},{"v"}}],FormatType->TraditionalForm]<>"=",ToString[MatrixForm[xMat],FormatType->StandardForm]<>ToString[MatrixForm[generalDegreesOfFreedom],FormatType->StandardForm]]

Now find the shape functions. Start by expression nodal unknowns in terms of nodal coordinates

(Subscript[u, 1])=u/.{x->(Subscript[x, 1]),y->(Subscript[y, 1])}
(Subscript[v, 1])=v/.{x->(Subscript[x, 1]),y->(Subscript[y, 1])}
(Subscript[u, 2])=u/.{x->(Subscript[x, 2]),y->(Subscript[y, 2])}
(Subscript[v, 2])=v/.{x->(Subscript[x, 2]),y->(Subscript[y, 2])}
(Subscript[u, 3])=u/.{x->(Subscript[x, 3]),y->(Subscript[y, 3])}
(Subscript[v, 3])=v/.{x->(Subscript[x, 3]),y->(Subscript[y, 3])}

Write the u=A a equation

{b,aMat}=Normal[CoefficientArrays[{(Subscript[u, 1]),(Subscript[v, 1]),(Subscript[u, 2]),(Subscript[v, 2]),(Subscript[u, 3]),(Subscript[v, 3])},generalDegreesOfFreedom]];

Print[ToString[MatrixForm[Transpose[{{"(Subscript[u, 1])","(Subscript[v, 1])","(Subscript[u, 2])","(Subscript[v, 2])","(Subscript[u, 3])","(Subscript[v, 3])"}}]],FormatType->TraditionalForm]<>"=",ToString[MatrixForm[aMat],FormatType->StandardForm]<>ToString[MatrixForm[generalDegreesOfFreedom],FormatType->StandardForm]]

Find a=(A^(-1))u from the above by matrix inversion

shapeFunctions=xMat.Inverse[aMat];

Now find the B matrix from the above N matrix by multiplying the N matrix by the following differetial operartors matrix

oper={{D[#1,x]&,0&},{0&,D[#1,y]&},{D[#1,y]&,D[#1,x]&}};
Print[ToString[MatrixForm[oper],FormatType->TraditionalForm]];

Now find B = oper * N

bMat=Inner[#1[#2]&,oper,shapeFunctions,Plus];

(bMat=Simplify[Assuming[Element[{(Subscript[y, 1]),(Subscript[y, 2]),(Subscript[y, 3]),(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3])},Reals],bMat]])//MatrixForm

Factor the determinant term from the above to the outside.

den=Denominator[bMat[[1,1]]] ;
bMat=bMat*den;
Print[ToString[1/den,FormatType->StandardForm]<>ToString[MatrixForm[Simplify[bMat]],FormatType->StandardForm]]

But area of triangle is

area=(1/2) Cross[{(Subscript[x, 2])-(Subscript[x, 1]),(Subscript[y, 2])-(Subscript[y, 1]),0},{(Subscript[x, 3])-(Subscript[x, 1]),(Subscript[y, 3])-(Subscript[y, 1]),0}][[3]]

Hence B matrix becomes

Panel[Style[ToString[1/"2 area",FormatType->StandardForm]<>ToString[MatrixForm[Simplify[bMat]],FormatType->StandardForm],18]]

finalB=(1/(2 area))bMat;