An object modeled as a mesh. |

The continuity of the deformations across cell borders depends on the choice of the shape functions. The minimum required for physical plausibility is $G^0$, obtained with affine functions, as illustrated in the next figure, where sharp discontinuity is created in the middle of an initially straight line. Smoother continuity is possible using higher-degree functions.

FE continuity. Linear shape
functions create non-smooth deformations. |

The interpolation functions are constant over time. The state (position $x$, velocity $v$) of the simulated object is represented by the values at the control nodes:

$$

\begin{eqnarray*}

\ve{x} &=& (x_0, \cdots x_{n-1})^T \\

\ve{v} &=& (v_0, \cdots v_{n-1})^T

\end{eqnarray*}

$$

We use $X$ to denote the initial values, $u(t)$ to denote the displacements over time, and $$x(t)= X + u(t)$$ to denote the current values, as illustrated in the following figure.

At each point of the domain, the simulated values (position, velocity) are entirely defined by the values at the control nodes and the interpolation functions. Let $p$ be the position of a point:

$$ p(t) = \sum_{i=0}^2 w_i x_i(t) $$

where the interpolation values $w_i$ are computed based on the positions in the initial configuration.

A Finite Element. Left: the reference configuration. Right: the displaced configuration. |

The four displacement modes of a
material element. |

The local displacement of the material can be represented by an initially orthonormal basis, as shown in the previous figure, called the deformation gradient and denoted $F$. It corresponds to the relation between the current and initial coordinates:

$$ F = \pd{x}{X}$$

The deformation gradient is a matrix where each column represents one of the basis vectors. It is easy to compute in triangular elements, based on the edges. Let $\bar E = (\bar E_1 \bar E_2)$ be a matrix where each column vector is an edge $\bar E_i = X_i - X_0$ of the initial triangle. This undeformed state corresponds to the identity matrix $\bar E^{-1} \bar E $. During the simulation, we compute the basis matrix $E(t)$ which columns are the current edges $E_i(t) = x_i(t)-x_0(t)$, and the new basis is computed by replacing the initial $\bar E$ with the current edge matrix:

$$ F(t) = \bar E^{-1}E(t) $$

The deformation gradient is used to compute the strain, which is a measure of the deformation. Different strains have been defined, for different purposes. One of the most popular is based on the idea that the displacement is rigid if and only if the deformation gradient is an orthogonal matrix, i.e. $F^TF = I$, where $I$ is the identity matrix. The Green-Lagrange strain tensor is therefore defined as: $$ \strain = \frac{1}{2}( F^TF-I)$$

The strain tensor is a symmetric matrix, it has thus only 3 independent entries (or 6 in 3d). The Voigt notation puts these independent values in a vector:

$$

\hat \strain = \left( \begin{array}{ccc} \strain_{xx} & \strain_{yy} & \strain_{xy} \end{array} \right)^T

$$

$$ \mathcal{w} = \frac{1}{2} \hat \strain^T H \hat \strain

$$

where matrix $H$ encodes the material properties:

$$

H = \frac{E}{1-\nu^2} \left( \begin{array}{ccc}

1 & \nu & 0 \\

\nu & 1 & 0 \\

0 & 0 & 1-\nu

\end{array} \right)

$$

where $E$ is called Young's modulus and encodes the stiffness of the material (the higher the more difficult to deform), and $\nu$ is called Poisson's ratio and encodes the material's incompressible (between 0: compressible and 1:incompressible).

The energy density is a constant across the element, thus the total deformation energy is the volume $\vol$ of the undeformed element times the energy density:

$$

\energy = \vol \energydensity

$$

and the elastic force applied to the control nodes is, as usual, derived from the energy:

$$

\ve f^T = -\pd{\energy}{\ve x} = -\pd{\energy}{\strain} \pd{\strain}{\ve x} = -\vstrain H \pd{\vstrain}{\ve x}

$$

where vector $\ve f$ concatenates the forces applied to the vertices and $\ve x$ concatenates their positions. The forces applied by the element to its nodes are split and added to the global force vector based on the node indices in the mesh, as illustrated in the following figure.

Force vector assembly. The force
vector of an element (dark blue) is split and accumulated to the global
force vector (light blue) based on the node indices. |

The stiffness matrix is the derivative of the force with respect to the positions:

$$ \stma = \pd{\ve f}{\ve x} $$

For a triangle element in 2d, it is a $6\times 6$ symmetric matrix. For a whole mesh composed of $n$ vertices, it is a $2n \times 2n$ symmetric matrix computed by assembling all the element matrices. Each node of the element has an index in the global mesh. The element matrix is split and accumulated in the global matrix, as illustrated in the following figure.

Stiffness matrix assembly. The
stiffness matrix of an element (dark blue) is split and accumulated to
the global stiffness matrix (light blue) based on the node indices. |

Francois Faure, University of Grenoble. Main page