# Is it possible to give the closed-form of the stiffness matrix of triangular prism element?

Suppose the vertices of triangular prism are 1, 2, 3 (bottom triangle), 4, 5, 6 (top triangle), the coordinates of the vertices are $$(x_i, y_i, z_i), i=1,2,3,4,5,6.$$

We can write the stiffness matrix of FEM of the triangular prism as follows:

$$K_{i,j} = \int_0^1\int_0^{1-\xi}\int_{-1}^1 \begin{pmatrix}\frac{\eth N_i}{\eth\xi} & \frac{\eth N_i}{\eth\eta} & \frac{\eth N_i}{\eth\zeta}\end{pmatrix} (\textbf{J}^{-1})^T\textbf{J}^{-1}\begin{pmatrix}\frac{\eth N_j}{\eth\xi} \\ \frac{\eth N_j}{\eth\eta} \\ \frac{\eth N_j}{\eth\zeta}\end{pmatrix}det(\mathbf{J})d\xi d\eta d\zeta\,\\ i,j=1,2,...,6.$$

where

$$\textbf{J}=\begin{pmatrix}\frac{\eth x}{\eth\xi} & \frac{\eth y}{\eth\xi} & \frac{\eth z}{\eth\xi} \\ \frac{\eth x}{\eth\eta} & \frac{\eth y}{\eth\eta} & \frac{\eth z}{\eth\eta} \\ \frac{\eth x}{\eth\zeta} & \frac{\eth y}{\eth\zeta} & \frac{\eth z}{\eth\zeta} \end{pmatrix}$$

and $$\textbf{N}=\begin{pmatrix}N_1 & N_2 & N_3 & N_4 & N_5 & N_6\end{pmatrix}^T\\ \left\{\begin{array}{} N_1=\frac{1}{2}(1-\zeta)(1-\xi-\eta)\\ N_2=\frac{1}{2}\xi(1-\zeta)\\ N_3=\frac{1}{‌​2}\eta(1-\zeta)\\ N_4=\frac{1}{2}(1+\zeta)(1-\xi-‌​‌​\eta)\\ N_5=\frac{1}{2}\xi(1+\zeta)\\ N_6=\frac{1}{2}\eta(1+\zeta) \end{array} \right.\ \\ \left\{\begin{array}{} x=\mathbf{N^T}\mathbf{x_e}\\ y=\mathbf{N^T}\mathbf{y_e}\\ z=\mathbf{N^T}\mathbf{z_e}\\ \end{array} \right.\ \\ \left\{\begin{array}{} \mathbf{x_e}=\begin{pmatrix}x_1 & x_2 & x_3 & x_4 & x_5 & x_6\end{pmatrix}^T\\ \mathbf{y_e}=\begin{pmatrix}y_1 & y_2 & y_3 & y_4 & y_5 & y_6\end{pmatrix}^T\\ \mathbf{z_e}=\begin{pmatrix}z_1 & z_2 & z_3 & z_4 & z_5 & z_6\end{pmatrix}^T \end{array} \right.\$$

Is it possible to give the closed-form of the integration by the powerful mathematica?

Thanks, Tang Laoya

-
What does $\zeta_i$ stand for? You have $\zeta$ as a scalar integration variable, so it can't be referring to vector components. – DumpsterDoofus Feb 9 '14 at 23:58
oh, sorry, $$\zeta_i$$ are the z-coordinate of standard triangular prism, i.e., $$\zeta_i=-1, i=1,2,3, \\ \zeta_i=1, i=4,5,6.$$ Thanks – Tang Laoya Feb 10 '14 at 0:07
Ok, so $\zeta_i$ is the $z$-coordinate of the $i$th vertex of the prism? Are the $x$ and $y$-coordinates of the vertices not used anywhere, only the $z$-components are used? And what are the $x_i,y_i,z_i$? Are they arbitrary constants? And what does $J^{-T}$ represent? Is it shorthand for $(J^{-1})^T$? – DumpsterDoofus Feb 10 '14 at 0:55
I have a feeling it won't be possible to get an analytic form, since $J^{-1}$ is going to have entries which are fractions that have polynomials in their denominators, which might be problematic. But I might be wrong. – DumpsterDoofus Feb 10 '14 at 1:01
Hi DumpsterDoofus, thanks for your kindly reply. The $$x_i, y_i, z_i$$ are coordinates of vertices of original triangular prism. The x and y-coordinates of the vertices are used in the integration. Yes, $$J^{-T}$$ is the shorthand for $$(J^{-1})^{T}$$Thanks – Tang Laoya Feb 10 '14 at 1:02

So with the following code:

n = 1/2 {(1 - \[Zeta]) (1 - \[Xi] - \[Eta]), \[Xi] (1 - \[Zeta]), \
\[Eta] (1 - \[Zeta]), (1 + \[Zeta]) (1 - \[Xi] - \[Eta]), \[Xi] (1 + \
\[Zeta]), \[Eta] (1 + \[Zeta])};
xe = {x1, x2, x3, x4, x5, x6};
ye = {y1, y2, y3, y4, y5, y6};
ze = {z1, z2, z3, z4, z5, z6};
x = n.xe;
y = n.ye;
z = n.ze;
J = D[{x, y, z}, {{\[Xi], \[Eta], \[Zeta]}, 1}]\[Transpose];
s[i_] := {D[n[[i]], \[Xi]], D[n[[i]], \[Eta]], D[n[[i]], \[Zeta]]};
k[i_, j_] := s[i].Inverse[J.J\[Transpose]].s[j] Det[J];


Executing k[1,1] yields an enormous expression which is 38 pages long when printed. Integrating it would be extremely difficult, unless there is some hidden way to simplify it. So I'm not sure it's possible.

-
Thanks for your so beautiful code:) – Tang Laoya Feb 10 '14 at 2:28
Could you please also give me the code to output the whole matrix K by integrating? Thanks – Tang Laoya Feb 10 '14 at 2:34
Yeah, I wish it gave a useful answer but unfortunately it doesn't. Kind of surprising; did you say that the tetrahedron matrix had a closed-form? I'd be curious to see how they calculated it, and why a tetrahedron is so much less complicated than a triangular prism. – DumpsterDoofus Feb 10 '14 at 2:41
Thanks for your kindly reply. There are many papers about tetrahedron, for example: dspace.uta.edu/bitstream/handle/10106/933/… – Tang Laoya Feb 10 '14 at 2:48
Well, if the tetrahedron case took an entire PhD thesis to find, I imagine the triangular prism case will be similarly hard. – DumpsterDoofus Feb 10 '14 at 4:39