sofa::component::forcefield::TetrahedronFEMForceFieldInternalData< gpu::cuda::CudaVectorTypes< TCoord, TDeriv, TReal > >::GPUElement Struct Reference

Static data associated with each element. More...

#include <CudaTetrahedronFEMForceField.h>

Collaboration diagram for sofa::component::forcefield::TetrahedronFEMForceFieldInternalData< gpu::cuda::CudaVectorTypes< TCoord, TDeriv, TReal > >::GPUElement:

[legend]


Data Fields
index of the 4 connected vertices
int	ia [BSIZE]
int	ib [BSIZE]
int	ic [BSIZE]
int	id [BSIZE]
material stiffness matrix
Real	gamma_bx2 [BSIZE]
Real	mu2_bx2 [BSIZE]
initial position of the vertices in the local (rotated) coordinate system
Real	bx [BSIZE]
Real	cx [BSIZE]
Real	cy [BSIZE]
Real	dx [BSIZE]
Real	dy [BSIZE]
Real	dz [BSIZE]

Real	Jbx_bx [BSIZE]
	strain-displacement matrix
Real	Jby_bx [BSIZE]
Real	Jbz_bx [BSIZE]

Detailed Description

template<class TCoord, class TDeriv, class TReal>
struct sofa::component::forcefield::TetrahedronFEMForceFieldInternalData< gpu::cuda::CudaVectorTypes< TCoord, TDeriv, TReal > >::GPUElement

Static data associated with each element.

Only supports the LARGE computation method, allowing for several simplifications.

For each element ABCD, the original algorithm computes the $3 \times 3$ rotation matrix R as follow:

the X axis is aligned to the AB edge
the Z axis is aligned to the cross-product between AB and AC
the Y axis is aligned to the cross-product between the previous two edges

This matrix is used to compute the initial position of each element in the local rotated coordinate frame, with the origin on A, as follow:

$\begin{eqnarray*} a &=& R^t \cdot (A-A) = \left( 0 \quad 0 \quad 0 \right) \\ b &=& R^t \cdot (B-A) = \left( b_x \quad 0 \quad 0 \right) \\ c &=& R^t \cdot (C-A) = \left( c_x \quad c_y \quad 0 \right) \\ d &=& R^t \cdot (D-A) = \left( d_x \quad d_y \quad d_z \right) \\ \end{eqnarray*}$

The material stiffness is handled by a $6 \times 6$ matrix K computed as follow:

$\begin{eqnarray*} K &=& \begin{pmatrix} \gamma+\mu^2 & \gamma & \gamma & 0 & 0 & 0 \\ \gamma & \gamma+\mu^2 & \gamma & 0 & 0 & 0 \\ \gamma & \gamma & \gamma+\mu^2 & 0 & 0 & 0 \\ 0 & 0 & 0 & \mu^2/2 & 0 & 0 \\ 0 & 0 & 0 & 0 & \mu^2/2 & 0 \\ 0 & 0 & 0 & 0 & 0 & \mu^2/2 \\ \end{pmatrix} \\ \gamma &=& \frac{youngModulus}{6 vol(ABCD)}\frac{poissonRatio}{(1 + poissonRatio)(1 - 2 poissonRatio)} \\ \mu^2 &=& \frac{youngModulus}{6 vol(ABCD)}\frac{1}{1+poissonRatio} \end{eqnarray*}$

The $12 \times 6$ Strain-displacement matrix is computed from the rotated initial positions abcd as follow:

$\begin{eqnarray*} J &=& \begin{pmatrix} -p_x(bcd) & 0 & 0 & -p_y(bcd) & 0 & -p_z(bcd) \\ 0 & -p_y(bcd) & 0 & -p_x(bcd) & -p_z(bcd) & 0 \\ 0 & 0 & -p_z(bcd) & 0 & -p_y(bcd) & -p_x(bcd) \\ p_x(cda) & 0 & 0 & p_y(cda) & 0 & p_z(cda) \\ 0 & p_y(cda) & 0 & p_x(cda) & p_z(cda) & 0 \\ 0 & 0 & p_z(cda) & 0 & p_y(cda) & p_x(cda) \\ -p_x(dab) & 0 & 0 & -p_y(dab) & 0 & -p_z(dab) \\ 0 & -p_y(dab) & 0 & -p_x(dab) & -p_z(dab) & 0 \\ 0 & 0 & -p_z(dab) & 0 & -p_y(dab) & -p_x(abc) \\ p_x(abc) & 0 & 0 & p_y(abc) & 0 & p_z(abc) \\ 0 & p_y(abc) & 0 & p_x(abc) & p_z(abc) & 0 \\ 0 & 0 & p_z(abc) & 0 & p_y(abc) & p_x(abc) \\ \end{pmatrix} \\ p(uvw) &=& u \times v + v \times w + w \times u \end{eqnarray*}$

Given the zeros in abcd, we have:

$\begin{eqnarray*} Jb &=& p(cda) = c \times d + d \times a + a \times c = \begin{pmatrix}c_x \\ c_y \\ 0\end{pmatrix} \times \begin{pmatrix}d_x \\ d_y \\ d_z\end{pmatrix} = \begin{pmatrix}c_y d_z \\ - c_x d_z \\ c_x d_y - c_y d_x\end{pmatrix} \\ Jc &=& -p(dab) = - d \times a - a \times b - b \times d = - \begin{pmatrix}b_x \\ 0 \\ 0\end{pmatrix} \times \begin{pmatrix}d_x \\ d_y \\ d_z\end{pmatrix} = \begin{pmatrix}0 \\ b_x d_z \\ - b_x d_y\end{pmatrix} \\ Jd &=& p(abc) = a \times b + b \times c + c \times a = \begin{pmatrix}b_x \\ 0 \\ 0\end{pmatrix} \times \begin{pmatrix}c_x \\ c_y \\ 0\end{pmatrix} = \begin{pmatrix}0 \\ 0 \\ b_x c_y\end{pmatrix} \end{eqnarray*}$

Also, as the sum of applied forces must be zero, we know that:

$\begin{eqnarray*} Ja+Jb+Jc+Jd &=& 0 \\ -p(bcd)+p(cda)-p(dab)+b(abc) &=& - b \times c - c \times d - d \times b \\ &+& c \times d + d \times a + a \times c \\ &-& d \times a - a \times b - b \times d \\ &+& a \times b + b \times c + c \times a \\ &=& 0 \\ Ja &=& -Jb-Jc-Jd \end{eqnarray*}$

The forces will be computed using $F = R J K J^t R^t X$ . We can apply a scaling factor to J if we divide K by its square: $F = R (1/b_x J) (b_x^2 K) (1/b_x J)^t R^t X$ .

This allows us to do all computations from the values $\left(b_x \quad c_x \quad c_y \quad d_x \quad d_y \quad d_z \quad \gamma b_x^2 \quad \mu^2 b_x^2 \quad Jb_x/b_x \quad Jb_y/b_x \quad Jb_z/b_x\right)$ . Including the 4 vertex indices, this represents 15 values.