% _________________________________________________________________________
%| 20NODE-60DOF HEXAHEDRAL SOLID Shape function control [A.] (C)(R) |
%|########################################## |
%| MECHANISM | Edge Upperside 4X(1Node-3Dof)[ux,vy,wz] |
%| | Underside 4X(1Node-3Dof)[ux,vy,wz] |
%| | Middle node N/A |
%| | Centroidal node N/A |
%|---------------------|---------------------------------------------------|
%| ANALYSIS | Shape function:Parametric (Pascal triangle) |
%| | 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. |
%|_________________________________________________________________________|
% | MATLAB 7.1 |
clc
clear
syms x y z real
%Element displacement function polynomial terms
w=[1,x,y,z,x^2,y^2,z^2,x*y,y*z,x*z,x^2*y,x^2*z,y^2*x,y^2*z,z^2*x,z^2*y,x*y*z,x^2*y*z,y^2*x*z,z^2*x*y];
Bound=[
subs(w,{x,y,z},{-1,-1,-1})
subs(w,{x,y,z},{ 1,-1,-1})
subs(w,{x,y,z},{ 1, 1,-1})
subs(w,{x,y,z},{-1, 1,-1})
subs(w,{x,y,z},{ 0,-1,-1})
subs(w,{x,y,z},{ 1, 0,-1})
subs(w,{x,y,z},{ 0, 1,-1})
subs(w,{x,y,z},{-1, 0,-1})
subs(w,{x,y,z},{-1,-1,0})
subs(w,{x,y,z},{ 1,-1,0})
subs(w,{x,y,z},{ 1, 1,0})
subs(w,{x,y,z},{-1, 1,0})
subs(w,{x,y,z},{-1,-1,1})
subs(w,{x,y,z},{ 1,-1,1})
subs(w,{x,y,z},{ 1, 1,1})
subs(w,{x,y,z},{-1, 1,1})
subs(w,{x,y,z},{ 0,-1,1})
subs(w,{x,y,z},{ 1, 0,1})
subs(w,{x,y,z},{ 0, 1,1})
subs(w,{x,y,z},{-1, 0,1})];
%Element node freedom displacements
syms u1 u2 u3 u4 u5 u6 u7 u8 u9 u10 u11 u12 u13 u14 u15 u16 u17 u18 u19 u20 real
dep=[u1,u2,u3,u4,u5,u6,u7,u8,u9,u10,u11,u12,u13,u14,u15,u16,u17,u18,u19,u20];
fun=w*Bound^-1*dep';
%Element shape functions
Shape=[diff(fun,u1,1)
diff(fun,u2,1)
diff(fun,u3,1)
diff(fun,u4,1)
diff(fun,u5,1)
diff(fun,u6,1)
diff(fun,u7,1)
diff(fun,u8,1)
diff(fun,u9,1)
diff(fun,u10,1)
diff(fun,u11,1)
diff(fun,u12,1)
diff(fun,u13,1)
diff(fun,u14,1)
diff(fun,u15,1)
diff(fun,u16,1)
diff(fun,u17,1)
diff(fun,u18,1)
diff(fun,u19,1)
diff(fun,u20,1)];
%Know shape function collocated finding shape functions
a=2;b=2;c=2;
Shapecontrol=[1/8*(1-2*x/a)*(1-2*y/b)*(1-2*z/c)*(-2*x/a-2*y/b-2*z/c-2)
1/8*(1+2*x/a)*(1-2*y/b)*(1-2*z/c)*(+2*x/a-2*y/b-2*z/c-2)
1/8*(1+2*x/a)*(1+2*y/b)*(1-2*z/c)*(+2*x/a+2*y/b-2*z/c-2)
1/8*(1-2*x/a)*(1+2*y/b)*(1-2*z/c)*(-2*x/a+2*y/b-2*z/c-2)
1/4*(1-(2*x/a)^2)*(1-2*y/b)*(1-2*z/c)
1/4*(1-(2*y/b)^2)*(1+2*x/a)*(1-2*z/c)
1/4*(1-(2*x/a)^2)*(1+2*y/b)*(1-2*z/c)
1/4*(1-(2*y/b)^2)*(1-2*x/a)*(1-2*z/c)
1/4*(1-(2*z/c)^2)*(1-2*x/a)*(1-2*y/b)
1/4*(1-(2*z/c)^2)*(1+2*x/a)*(1-2*y/b)
1/4*(1-(2*z/c)^2)*(1+2*x/a)*(1+2*y/b)
1/4*(1-(2*z/c)^2)*(1-2*x/a)*(1+2*y/b)
1/8*(1-2*x/a)*(1-2*y/b)*(1+2*z/c)*(-2*x/a-2*y/b+2*z/c-2)
1/8*(1+2*x/a)*(1-2*y/b)*(1+2*z/c)*(+2*x/a-2*y/b+2*z/c-2)
1/8*(1+2*x/a)*(1+2*y/b)*(1+2*z/c)*(+2*x/a+2*y/b+2*z/c-2)
1/8*(1-2*x/a)*(1+2*y/b)*(1+2*z/c)*(-2*x/a+2*y/b+2*z/c-2)
1/4*(1-(2*x/a)^2)*(1-2*y/b)*(1+2*z/c)
1/4*(1-(2*y/b)^2)*(1+2*x/a)*(1+2*z/c)
1/4*(1-(2*x/a)^2)*(1+2*y/b)*(1+2*z/c)
1/4*(1-(2*y/b)^2)*(1-2*x/a)*(1+2*z/c)];
factor(Shape-Shapecontrol)
%IF all terms equal zero than [OK]
