Discretiztion

The weak form of principle of virtual work (PVW) is discretized by using a shape function to map the distribution of the displacement to the whole cell. The shape function is based on the Lagrangian polynominals. The Lagrangian polynominal interpolates the the displacement between discrete nodes located inside the cell also referred to as element. The number of nodes per element is defined by the p-refinement, which is explained in the next section. The minimum number of nodes is 2 in each space direction. Hence, in 3D, 8 nodes (one at each element vertex) exist per element. The Lagrangian polynominal reads

\begin{equation} N_{ijk}(x, y, z) = \prod_{0 \leq l \leq (p+2); l \neq i} \frac{x - x_i}{x_l - x_i} \prod_{0 \leq l \leq (p+2); l \neq j} \frac{y - y_j}{y_l - y_j} \prod_{0 \leq l \leq (p+2); l \neq k} \frac{z - z_k}{z_l - z_k} \end{equation}

The subscribts \( i\), \( j\), and \( k \) denote the node position in \( x \), \( y \), and \( z \) direction. That is, \( N_{0,0,0} \) denotes the zeroth node position in \( x \), \( y \), and \( z \) direction. \( N_{1,2,1} \) denotes the first node position in \( x \), the second node position in \( y \) and the first node position in \( z \) direction. The separte shape functions are assembled in the matrix \( \mv{\mv{N}} \). To obtain the strain-interpolation matrix \( \mv{\mv{B}} \), the derivative of matrix \( \mv{\mv{N}} \) is calculated. The linear system of equations solved by the FC solver reads

\begin{equation} \mv{\mv{K}} \mv{u} = \mv{F} \end{equation}

where \( \mv{\mv{K}} \) is the gobal stiffness matrix, \( \mv{u} \) is the global displacement vector, and \( \mv{F} \) is the global force vector. The global vectors and matrices are assembled from the element matrices/vectors. The element stiffness matrix is defined by

\begin{equation} \mv{\mv{K}} = \int\limits_{V_e} \mv{\mv{B}}^T \mv{\mv{C}} \mv{\mv{B}} det|\mv{\mv{X}}| \, dV_e \end{equation}

where \( V_e \) is the element volume and \( det|\mv{\mv{X}}| \) is the determinant of the Jacobian matrix.

P-refinement

When using p-refinement, the number of nodes per element is increased. That is, increasing the p-refinement by \( 1 \) increases the number of nodes per Cartesian direction by \( 1 \). The total number of nodes per element is given by \( (p + 2)(p + 2)(p + 2) = p^3 + 6p^2 + 12p + 8 \). The distribution of the points follows the distribution of the Lobatto-Points or an equidistant distribution.

Warning: Although the p-refinement is implemented until p = 10, it is tested so far only for upto p = 3.

Subcell integration

To increase the accuracy at boundaries without increasing the degree of freedom, a subcell integration can be applied at curved boundaries. That is, the boundary cell is subdivided into \( 4 \) (2D) or \( 8 \) (3D) cells. For each subcell it is checked, if the subcell has a cut with the geometry. If yes, the cell is further subdivided. If not, the integration is performed on the cell using \( \alpha \ll 1 \) if the cell is located outside and \( \alpha = 1 \) if the cell is inside the geometry. The subdivision is performed until the maximum subcell depth is reached. On the maximum depth, the inside-outside check is conducted for each integration point and the \( \alpha \) is changed depending on the check.

Note: Although the degree of freedome is not increased due to the subcell integration, the computational cost are increased. For example, a subcell depth of 10 results in upto \( 10^9 \) subcells per boundary cell.

Boundary conditions

In the following, the two types of boundaries are shown.

Essential boundary conditions

At essential boundaries the displacement is set. This is done by solving the following integral and adding the resulting matrix to the element stiffness matrix.

\begin{equation} \mv{\mv{K}} = \mv{\mv{K}} + \int\limits_A \mv{\mv{N}}^T k \mv{\mv{N}} det|\mv{\mv{X}}| \, dA \end{equation}

Here, \( k \) is a factor with \( k \gg 1 \) and \( A \) is the boundary surface, which is approximated by the triangles of the STL file. Essential boundary conditions have the following numbers, which can be set in the geometry.toml file:

BC.0 = 8010
 
BC.0 = 8011
 
BC.0 = 8012

Note: So far, only a displacement of \(u = 0\) can be prescribed.

Natural boundary conditions

At natural boundaries the surface traction is set. This is done by solving the following integral and adding the resulting vector to the element force vector.

\begin{equation} \mv{F} = \mv{F} + \int\limits_A \mv{\mv{N}}^T \mv{t} \, det|\mv{\mv{X}}| \, dA \end{equation}

Here, \( t \) denotes the surface traction and \( A \) is the boundary surface, which is approximated by the triangles of the STL file. Natural boundary conditions have the following numbers, which can be set in the geometry.toml file:

BC.0 = 8030
 
BC.0 = 8031
 
BC.0 = 8032
 
BC.0 = 3035