% _________________________________________________________________________
%|____10NODE_20DOF PARAMETRIC TRIANGULAR PLANE STRESS FEM_______[A.]_(C)(R)|
%| |
%| Shape function:Homogen |
%| Stifness :Topology+Accumulate method |
%|_______________________________Solve Equation:Cholesky-[L][D][u]_________|
%| Maxima 5.9.0.9beta2 http://maxima.sourceforge.net |
%| Using Lisp Kyoto Common Lisp GCL 2.6.3 (aka GCL) |
%| Distributed under the GNU Public License. See the file COPYING. |
%| Dedicated to the memory of William Schelter. |
%| This is a development version of Maxima. The function bug_report() |
%| provides bug reporting information. |
%| |
%|_____PARAMETRIC FEM ANALYSIS SUBPROGRMAM_______Maxima5.9.0.9beta2_[GPL]__|
%MAXIMA CODE
%========================================This procedure runing MAXIMA
%PROCEDURE(1):Information element shape function depend the cartesian
% coordinates for homogen axis.
%PROCEDURE(2):Element Stiffness matrix for Local element`s axis
%__________________________________________[PROCEDURE(1)];
%Alan:0.5*determinant(matrix([1,x1,y1],[1,x2,y2],[1,x3,y3]));
%fe1:determinant(matrix([1,x,y],[1,x2,y2],[1,x3,y3]))/(2*Alan);
%fe2:determinant(matrix([1,x,y],[1,x3,y3],[1,x1,y1]))/(2*Alan);
%fe3:determinant(matrix([1,x,y],[1,x1,y1],[1,x2,y2]))/(2*Alan);
%'e1=(-x*(Y3-Y2)+x2*Y3-x3*Y2+(x3-x2)*y)/(-x1*(Y3-Y2)+x2*Y3-x3*Y2+(x3-x2)*Y1);
%'e2=(-x1*Y3-x*(Y1-Y3)+x3*Y1+(x1-x3)*y)/(-x1*(Y3-Y2)+x2*Y3-x3*Y2+(x3-x2)*Y1);
%'e3=(-x(Y2-Y1)+x1*Y2-x2*Y1+(x2-x1)*y)/(-x1*(Y3-Y2)+x2*Y3-x3*Y2+(x3-x2)*Y1);
%'Area=0.5*(-x1*(Y3-Y2)+x2*Y3-x3*Y2+(x3-x2)*Y1);
%'Ai:( x2*Y3-x3*Y2);
%'Aj:(-x1*Y3+x3*Y1);
%'Ak:( x1*Y2-x2*Y1);
%'Bi:-(Y3-Y2);
%'Bj:-(Y1-Y3);
%'Bk:-(Y2-Y1);
%'Ci:+(x3-x2);
%'Cj:+(x1-x3);
%'Ck:+(x2-x1);
%___________________________________________[PROCEDURE(2)];
%fe1:1/(2*Area)*(Ai+Bi*x+Ci*y);
%fe2:1/(2*Area)*(Aj+Bj*x+Cj*y);
%fe3:1/(2*Area)*(Ak+Bk*x+Ck*y);
%Nh1(x,y):=1/2*e1*(3*e1-1)*(3*e1-2);
%Nh2(x,y):=1/2*e2*(3*e2-1)*(3*e2-2);
%Nh3(x,y):=1/2*e3*(3*e3-1)*(3*e3-2);
%Nh4(x,y):=9/2*e1*e2*(3*e1-1);
%Nh5(x,y):=9/2*e1*e2*(3*e2-1);
%Nh6(x,y):=9/2*e2*e3*(3*e2-1);
%Nh7(x,y):=9/2*e2*e3*(3*e3-1);
%Nh8(x,y):=9/2*e3*e1*(3*e3-1);
%Nh9(x,y):=9/2*e3*e1*(3*e1-1);
%Nh10(x,y):=27*e1*e2*e3;
%depends(f,[e1,e2,e3],[e1,e2,e3],[x,y]);
%diff(f,x);
%Hx1(x,y):=diff(Nh1(x,y),e1)*diff(fe1,x)+diff(Nh1(x,y),e2)*diff(fe2,x)+diff(Nh1(x,y),e3)*diff(fe3,x);
%Hx2(x,y):=diff(Nh2(x,y),e1)*diff(fe1,x)+diff(Nh2(x,y),e2)*diff(fe2,x)+diff(Nh2(x,y),e3)*diff(fe3,x);
%Hx3(x,y):=diff(Nh3(x,y),e1)*diff(fe1,x)+diff(Nh3(x,y),e2)*diff(fe2,x)+diff(Nh3(x,y),e3)*diff(fe3,x);
%Hx4(x,y):=diff(Nh4(x,y),e1)*diff(fe1,x)+diff(Nh4(x,y),e2)*diff(fe2,x)+diff(Nh4(x,y),e3)*diff(fe3,x);
%Hx5(x,y):=diff(Nh5(x,y),e1)*diff(fe1,x)+diff(Nh5(x,y),e2)*diff(fe2,x)+diff(Nh5(x,y),e3)*diff(fe3,x);
%Hx6(x,y):=diff(Nh6(x,y),e1)*diff(fe1,x)+diff(Nh6(x,y),e2)*diff(fe2,x)+diff(Nh6(x,y),e3)*diff(fe3,x);
%Hx7(x,y):=diff(Nh7(x,y),e1)*diff(fe1,x)+diff(Nh7(x,y),e2)*diff(fe2,x)+diff(Nh7(x,y),e3)*diff(fe3,x);
%Hx8(x,y):=diff(Nh8(x,y),e1)*diff(fe1,x)+diff(Nh8(x,y),e2)*diff(fe2,x)+diff(Nh8(x,y),e3)*diff(fe3,x);
%Hx9(x,y):=diff(Nh9(x,y),e1)*diff(fe1,x)+diff(Nh9(x,y),e2)*diff(fe2,x)+diff(Nh9(x,y),e3)*diff(fe3,x);
%Hx10(x,y):=diff(Nh10(x,y),e1)*diff(fe1,x)+diff(Nh10(x,y),e2)*diff(fe2,x)+diff(Nh10(x,y),e3)*diff(fe3,x);
%Hx1(x,y);
%Hx2(x,y);
%Hx3(x,y);
%Hx4(x,y);
%Hx5(x,y);
%Hx6(x,y);
%Hx7(x,y);
%Hx8(x,y);
%Hx9(x,y);
%Hx10(x,y);
%Hy1(x,y):=diff(Nh1(x,y),e1)*diff(fe1,y)+diff(Nh1(x,y),e2)*diff(fe2,y)+diff(Nh1(x,y),e3)*diff(fe3,y);
%Hy2(x,y):=diff(Nh2(x,y),e1)*diff(fe1,y)+diff(Nh2(x,y),e2)*diff(fe2,y)+diff(Nh2(x,y),e3)*diff(fe3,y);
%Hy3(x,y):=diff(Nh3(x,y),e1)*diff(fe1,y)+diff(Nh3(x,y),e2)*diff(fe2,y)+diff(Nh3(x,y),e3)*diff(fe3,y);
%Hy4(x,y):=diff(Nh4(x,y),e1)*diff(fe1,y)+diff(Nh4(x,y),e2)*diff(fe2,y)+diff(Nh4(x,y),e3)*diff(fe3,y);
%Hy5(x,y):=diff(Nh5(x,y),e1)*diff(fe1,y)+diff(Nh5(x,y),e2)*diff(fe2,y)+diff(Nh5(x,y),e3)*diff(fe3,y);
%Hy6(x,y):=diff(Nh6(x,y),e1)*diff(fe1,y)+diff(Nh6(x,y),e2)*diff(fe2,y)+diff(Nh6(x,y),e3)*diff(fe3,y);
%Hy7(x,y):=diff(Nh7(x,y),e1)*diff(fe1,y)+diff(Nh7(x,y),e2)*diff(fe2,y)+diff(Nh7(x,y),e3)*diff(fe3,y);
%Hy8(x,y):=diff(Nh8(x,y),e1)*diff(fe1,y)+diff(Nh8(x,y),e2)*diff(fe2,y)+diff(Nh8(x,y),e3)*diff(fe3,y);
%Hy9(x,y):=diff(Nh9(x,y),e1)*diff(fe1,y)+diff(Nh9(x,y),e2)*diff(fe2,y)+diff(Nh9(x,y),e3)*diff(fe3,y);
%Hy10(x,y):=diff(Nh10(x,y),e1)*diff(fe1,y)+diff(Nh10(x,y),e2)*diff(fe2,y)+diff(Nh10(x,y),e3)*diff(fe3,y);
%Hy1(x,y);
%Hy2(x,y);
%Hy3(x,y);
%Hy4(x,y);
%Hy5(x,y);
%Hy6(x,y);
%Hy7(x,y);
%Hy8(x,y);
%Hy9(x,y);
%Hy10(x,y);
%Af(x,y):=matrix([Hx1(x,y) , 0 ,Hx2(x,y) , 0 ,Hx3(x,y) , 0 ,Hx4(x,y) , 0 ,Hx5(x,y) , 0 ,Hx6(x,y) , 0 ,Hx7(x,y) , 0 ,Hx8(x,y) , 0 ,Hx9(x,y) , 0 ,Hx10(x,y) , 0 ],
% [ 0 ,Hy1(x,y) , 0 , Hy2(x,y) , 0 , Hy3(x,y) , 0 ,Hy4(x,y) , 0 , Hy5(x,y) , 0 ,Hy6(x,y), 0 , Hy7(x,y) , 0 ,Hy8(x,y) , 0 , Hy9(x,y) , 0 ,Hy10(x,y)],
% [Hy1(x,y) , Hx1(x,y) , Hy2(x,y) , Hx2(x,y) , Hy3(x,y) , Hx3(x,y) ,Hy4(x,y) ,Hx4(x,y) , Hy5(x,y) , Hx5(x,y) ,Hy6(x,y) ,Hx6(x,y) , Hy7(x,y) , Hx7(x,y) ,Hy8(x,y) ,Hx8(x,y) , Hy9(x,y) , Hx9(x,y) ,Hy10(x,y) ,Hx10(x,y)]);
%Af(x,y);
%A(x,y):=ratsubst((1-t)*(1-e1),e3,ratsubst((1-e1)*t,e2,Af(x,y)));
%A(x,y);
%C(x,y):=matrix([1 , v , 0 ],
% [v , 1 , 0 ],
% [0 , 0 ,(1-v)/2]);
%C(x,y);
%B(x,y):=Transpose(A(x,y));
%Ku(x,y):=B(x,y).C(x,y).A(x,y);
%K(x,y):=KSI*factor(defint(defint(2*Area*Ku(x,y)*(1-e1),e1,0,1),t,0,1));
%fortmx(K,K(x,y));
%===========Information:
%Element stiffness matrix high order triangular f.e.m for not
%integrated. This region requiremented Beta and Gamma partial functions
%this integration. Nevertless this integration not parametric terms and
%imtegration not need. But High order stiffness matrix terms it`s need.
%e1=e1
%e2=(1-e1)*t
%e3=(1-t)*(1-e1)
%===============================================END MAXIMA
