Tet4 element

Tet4 is the simplest linear-interpolation tetrahedron-shaped element with 4 nodes, each having only 3 translational degrees of freedom (rotations of nodes are ignored).

Element coordinate system

Initial element configuration is the reference one, and element-local nodal positions (orientations are ignored) are stored in refPos. This reference configuration is created automatically when the element is first used (usually when woo.dem.DynDt asks for stiffnesses before the simulation starts, or by an explicit call to setRefConf).

Current element configuration is stored in woo.fem.Tet4.node; its origin is always at element’s centorid (getCentroid), which is computed as average of nodal positions. The orientation is determined using the algorithm described in [MS14] (pg. 114, Procedure 1):

1. compute

   • \(C\), element-local reference coordinates of nodes (transpose of refPos);

   • \(D\), current globally-oriented coordinates of nodes, with the origin moved to the current element centroid (getCentroid)

2. we introduce the following notation:

   • non-unit quaternion \(\mathbf{q}\) is a 4-tuple \(\mathbf{q}=(q_0,q_x,q_y,q_z)=(q_0,q)\), where \(q_0\) is the real component and \((q_x,q_y,q_z)=q\) is 3D vector of imaginary components;

   • \(\mathbf{q}^\circ=(0,\hat q)\) forms pure-imaginary quaternion from \(q\);

   • \(\bullet_L\) and \(\bullet_R\) are matrices representing respectively left- and right-multiplication by quaternion \(\bullet\) ([MS14], A.4-A.5):

     \[ \begin{align}\begin{aligned}\begin{split}[\mathbf{q}]_L=\begin{pmatrix} q_0 & -q^T \\ q & q_0 I+\hat q\end{pmatrix},\end{split}\\\begin{split}[\mathbf{q}]_R=\begin{pmatrix} q_0 & -q^T \\ q & q_0 I-\hat q\end{pmatrix},\end{split}\end{aligned}\end{align} \]

     which simplifies, for pure-imaginary arguments, to

     \[ \begin{align}\begin{aligned}\begin{split}[\mathbf{q}^\circ]_L=\begin{pmatrix} 0 & -q^T \\ q & +\hat q\end{pmatrix},\end{split}\\\begin{split}[\mathbf{q}^\circ]_R=\begin{pmatrix} 0 & -q^T \\ q & -\hat q\end{pmatrix};\end{split}\end{aligned}\end{align} \]

   • the hat operator \(\hat\bullet\) is the skew-symmetric matrix representing cross-product by \(\bullet\) ([MS14], A.1):

     \[\begin{split}\hat a=\begin{pmatrix} 0 & -a_3 & a_2 \\ a_3 & 0 & -a_1 \\ -a_2 & a_1 & 0 \end{pmatrix}.\end{split}\]

3. compute [MS14], (8b):

   \[\mathcal{B}(C,D)=\sum_{n=1}^{4}\beta(c_n,d_n)^T\beta(c_n,d_n)\]

   where \(c_n\), \(d_n\) are respectively \(n\)-th rows of \(C\) and \(D\); those are 4-tuples, or non-unit quaternions \(\mathbf{q}=(q_0,q_x,q_y,q_z)\) and the \(\beta\) function is defined as

   \[\beta(c,d)=\left[d^\circ\right]_L-\left[c^\circ\right]_R\]

4. Compute the smallest eigenvalue of the 4x4 \(\mathcal{B}(C,D)\) matrix; since the matrix is symmetric, the eigenvalue is real; the corresponding eigenvector is normalized and assigned to ori of the element frame.

Stiffness matrix

Todo

Mostly from [Fel99a], chapter 9.

Lumped mass and inertia

Suppose the tetrahedron is defined by vertices \(\vec{v}_1, \vec{v}_2, \vec{v}_3, \vec{v}_4\).

Centroid of the tetrahedron is the average \(\vec{c}_g=\frac{1}{4}\vec{v}_i\).

Volume of tetrahedron can be computed in various ways, wikipedia gives the two following:

\[\begin{split}V=\frac{1}{6}\begin{vmatrix}(\vec{v}_1-\vec{v}_4)^T \\ (\vec{v}_2-\vec{v}_4)^T \\ (\vec{v}_3-\vec{v}_4)^T\end{vmatrix}=\frac{(\vec{v}_1-\vec{v}_4)\cdot\left((\vec{v}_2-\vec{v}_4)\times(\vec{v}_3-\vec{v}_4)\right)}{6}\end{split}\]

where non-canonical vertex ordering yields negative volume value.

Inertia tensor of tetrahedron in term of its vertex coordinates (with respect to origin and global axes) is derived at this excellent page. We only summarize the most important parts here.

In general, inertia tensor \(\mat{J}\) of any body can be computed from covariance

\[\mat{J}=\operatorname{tr}(\mat{C})\mat{I}_{3}-\mat{C}\]

where covariance is outer product of coordinates over the volume

\[\mat{C}\equiv\int_V \vec{x}\vec{x}^T \d V .\]

The page referenced shows that covariance of generic tetrahedron can be derived by transforming covariance of unit tetrahedron, giving

\[\mat{C}=\frac{V}{20}\left(\sum_i\vec{v}_i\vec{v}_i^T+\sum_i\vec{v}_i\sum\vec{v}_i^T\right) .\]

When lumping mass and inertia, only the part adjacent to each node is considered; tetrahedron is partitioned into 8 sub-tetrahedra. When we define \(\vec{v}_{ij}=\frac{\vec{v}_i+\vec{v}_j}{2}\), those sub-tetrahedra for node \(i\) are defined by vertices \(\{\vec{v}_i,\vec{v}_{ij},\vec{v}_{ik},\vec{v}_{il}\}\) and \(\{\vec{v}_{ij},\vec{v}_{ik},\vec{v}_{il},\vec{c}_g\}\) (partitioning planes are voronoi tesselation of vertices). Mass and inertia of sub-tetrahedra are computed in node coordinate system, summed over all attached particles, and lumped into the node.

Tip

Report issues or inclarities to github.