Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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

share|improve this question
    
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 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 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 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 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 at 1:02

1 Answer 1

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.

share|improve this answer
    
Thanks for your so beautiful code:) –  Tang Laoya Feb 10 at 2:28
    
Could you please also give me the code to output the whole matrix K by integrating? Thanks –  Tang Laoya Feb 10 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 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 at 2:48
1  
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 at 4:39

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.