6 Finite Element Models

6.1 Overview

The finite element method (FEM) is a numerical technique used for solving a system of partial differential equations (PDEs) over some domain. The general approach is to divide the domain into a set of building blocks, referred to as elements. These partition the space, and form local domains over which the system of equations can be locally approximated. The corners of these elements, the nodes, become control points in a discretized system. The solution is then assumed to be smoothly interpolated across the elements based on values determined at the nodes. Using this discretization, the differential system is converted into an algebraic one, which is often linearized and solved iteratively.

In ArtiSynth, the PDEs considered are the governing equations of continuum mechanics: the conservation of mass, momentum, and energy. To complete the system, a constitutive equation is required that describes the stress-strain response of the material. This constitutive equation is what distinguishes between material types. The domain is the three-dimensional space that the model occupies. This must be divided into small elements which accurately represent the geometry. Within each element, the PDEs are sampled at a set of points, referred to as integration points, and terms are numerically integrated to form an algebraic system to solve.

The purpose of the rest of this chapter is to describe the construction and use of finite elements models within ArtiSynth. It does not further discuss the mathematical framework or theory. For an in-depth coverage of the nonlinear finite element method, as applied to continuum mechanics, the reader is referred to the textbook by Bonet and Wood [3].

6.1.1 FemModel3d

The basic type of finite element model is implemented in the class FemModel3d. This class controls some properties that are used by the model as a whole. The key ones that affect simulation dynamics are:

Property Description
density The density of the model
material An object that describes the material’s constitutive law (i.e. its stress-strain relationship).
particleDamping Proportional damping associated with the particle-like motion of the FEM nodes.
stiffnessDamping Proportional damping associated with the system’s stiffness term.

These properties can be set and retrieved using the methods

    setDensity (double density);    // sets the density
    double getDensity ();           // gets the density
    setMaterial (FemMaterial mat);  // sets the FEM’s material
    FemMaterial getMaterial ();     // gets the FEM’s material
    setParticleDamping (double d);  // sets the particle (mass) damping coefficient
    double getParticleDamping ();   // gets the particle (mass) damping coefficient
    setStiffnessDamping (double d); // sets the stiffness damping coefficient
    double getStiffnessDamping ();  // gets the stiffness damping coefficient

Keep in mind that ArtiSynth is essentially “unitless” (Section 4.2), so it is the responsibility of the developer to ensure that all properties are specified in a compatible way.

The density of the model is used to compute the mass distribution throughout the volume. Note that we use a diagonally lumped mass matrix (DLMM) formulation, so the mass is assumed to be concentrated at the location of the discretized FEM nodes. To allow for a spatially-varying density, densities can be explicitly set for individual elements, or masses can be explicitly set for indiviual nodes.

The FEM’s material property is a delegate object used to compute stress and stiffness within individual elements. It handles the constitutive component of the model, as described in more detail in Section 6.1.3. In addition to the main material defined for the model, it is also possible set a material on a per-element basis, and to define additional materials which augment the behavior of the main materials (Section 6.8).

The two damping parameters are related to Rayleigh damping, which is used to dissipate energy within finite element models. There are two proportional damping terms: one related to the system’s mass, and one related to stiffness. The resulting damping force applied is

\displaystyle{\bf f}_{d} \displaystyle=-(d_{M}{\bf M}+d_{K}{\bf K}){\bf v}, (6.1)

where d_{M} is the value of particleDamping, d_{K} is the value of stiffnessDamping, {\bf M} is the FEM model’s lumped mass matrix, {\bf K} is the FEM’s stiffness matrix, and {\bf v} is the concatenated vector of FEM node velocities. Since the lumped mass matrix is diagonal, the mass-related component of damping can be applied separately to each FEM node. Thus, the mass component reduces to the same system as Equation (3.3), which is why it is referred to as “particle damping”.

6.1.2 Component Structure

Each FemModel3d contains several lists of sub-components:

nodes

The particle-like dynamic components of the model. These lie at the corners of the elements and carry all the mass (due to DLMM formulation).

elements

The volumetric model elements. These define the 3D sub-units over which the system is numerically integrated.

shellElements

The shell elements. These define additional 2D sub-units over which the system is numerically integrated.

meshes

The geometry in the model. This includes the surface mesh, and any other embedded geometries.

materials

Optional additional materials which can be added to the model to augment the behavior of the model’s material property. This is described in more detail in Section 6.8.

fields

Optional field components which can be used to interpolate application-defined quantities over the FEM model’s domain. Fields are described in detail in Section 6.10.

The nodes, elements and meshes components are illustrated in Figure 6.2.

(a) FEM model (b) Nodes (c) Elements (d) Geometry
Figure 6.2: Sub-components of FemModel3d.

6.1.2.1 Nodes


The set of nodes belong to a finite element model can be obtained by the method

PointList<FemNode3d> getNodes();  // returns list of FEM nodes

Nodes are implemented in the class FemNode3d, which is a subclass of Particle (Section 3.1). They are the main dynamic components of the finite element model. The key properties affecting simulation dynamics are:

Property Description
restPosition The initial position of the node.
position The current position of the node.
velocity The current velocity of the node.
mass The mass of the node.
dynamic Whether the node is considered dynamic or parametric (e.g. boundary condition).

Each of these properties has corresponding getXxx() and setXxx(...) functions to access and modify them.

The restPosition property defines the node’s position in the FEM model’s “natural” or “undeformed” state. Rest positions are used to compute an initial configuration for the model, from which strains are determined. A node’s rest position can be updated in code using the method: FemNode3d.setRestPosition(Point3d).

If any node’s rest positions are changed, the current values for stress and stiffness will become invalid. They can be manually updated using the method FemModel3d.updateStressAndStiffness() for the parent model. Otherwise, stress and stiffness will be automatically updated at the beginning of the next time step.

The properties position and velocity give the node’s current 3D state. These are common to all point-like particles, as is the mass property. Here, however, mass represents the lumped mass of the immediately surrounding material. Its value is initialized by equally dividing mass contributions from each adjacent element, given their densities. For a finer control of spatially-varying density, node masses can be set manually after FEM creation.

The FEM node’s dynamic property specifies whether or not the node is considered when computing the dynamics of the system. If not, it is treated as being parametrically controlled. This has implications when setting boundary conditions (Section 6.1.4).

6.1.2.2 Elements


Elements are the 3D volumetric spatial building blocks of the domain. Within each element, the displacement (or strain) field is interpolated from displacements at nodes:

\displaystyle{\bf u}({\bf x}) \displaystyle=\sum_{i=1}^{N}\phi_{i}({\bf x}){\bf u}_{i}, (6.2)

where {\bf u}_{i} is the displacement of the ith node that is adjacent to the element, and \phi_{i}(\cdot) is referred to as the shape function (or basis function) associated with that node. Elements are classified by their shape, number of nodes, and shape function order (Table 6.1). ArtiSynth supports the following element types:

TetElement, PyramidElement, WedgeElement, HexElement,
QuadtetElement QuadpyramidElement QuadwedgeElement QuadhexElement
Table 6.1: Supported element types
Element Type # Nodes Order # Integration Points
TetElement 4 linear 1
PyramidElement 5 linear 5
WedgeElement 6 linear 6
HexElement 8 linear 8
QuadtetElement 10 quadratic 4
QuadpyramidElement 13 quadratic 5
QuadwedgeElement 15 quadratic 9
QuadhexElement 20 quadratic 14

The base class for all of these is FemElement3d. A numerical integration is performed within each element to create the (tangent) stiffness matrix. This integration is performed by evaluating the stress and stiffness at a set of integration points within each element, and applying numerical quadrature. The list of elements in a model can be obtained with the method

  RenderableComponentList<FemElement3d> getElements();

All objects of type FemElement3d have the following properties:

Property Description
density Density of the element
material An object that describes the constitutive law within the element (i.e. its stress-strain relationship).

If left unspecified, the element’s density is inherited from the containing FemModel3d object. When set, the mass of the element is computed and divided amongst all its nodes, updating the lumped mass matrix.

Each element’s’ material property is also inherited by default from the containing FemModel3d. Specifying a material here allows for spatially-varying material properties across the model. Materials will be discussed further in Section 6.1.3.

6.1.2.3 Shell elements

Shell elements are additional 2D spatial building blocks which can be added to a model. They are typically used to model structures which are too thin to be easily represented by 3D volumetric elements, or to provide additional internal stiffness within a set of volumetric elements.

ArtiSynth presently supports the following shell element types, with the number of nodes, shape function order, and integration point count described in Table 6.2:

Table 6.2: Supported shell element types
Element Type # Nodes Order # Integration Points
ShellTriElement 3 linear 9 (3 if membrane)
ShellQuadElement 4 linear 8 (4 if membrane)

The base class for all shell elements is ShellElement3d, which contains the same density and material properties as FemElement3d, as well as the additional property defaultThickness, whose use will be described below.

The list of shell elements in a model can be obtained with the method

  RenderableComponentList<ShellElement3d> getShellElements();

Both the volumetric elements (FemElement3d) and the shell elements (ShellElement3d) derive from the base class FemElement3dBase. To obtain all the elements in an FEM model, both shell and volumetric, one may use the method

  ArrayList<FemElement3dBase> getAllElements();

Each shell element can actually be instantiated in two forms:

  • As a regular shell element, which has a bending stiffness;

  • As a membrane element, which does not have bending stiffness.

Regular shell elements are implemented using the same extensible director formulation used by FEBio [9], and more specifically the front/back node formulation [10]. Each node associated with a (regular) shell element is assigned a director, which is a 3D vector providing a normal direction and virtual thickness at that node (Figure 6.3). This virtual thickness allows us to continue to use 3D materials to provide the constitutive laws that determine the shell’s stress/strain response, including its bending behavior. It also allows us to continue to use the element’s density to determine its mass.

Figure 6.3: ShellQuadElement with the directors (dark blue lines) visible at the nodes.

Director information is automatically assigned to a FemNode3d whenever one or more regular shell elements is connected to it. This information includes both the current value of the director, its rest value, and its velocity, with the difference between the first two determining the element’s bending strain. These quantities can be queried using the methods

  Vector3d getDirector ();        // return the current director value
  void getDirector (Vector3d dir);
  Vector3d getRestDirector ();    // return the rest director value
  void getRestDirector (Vector3d dir);
  Vector3d getDirectorVel ();     // return the director velocity
  boolean hasDirector();          // does this node have a director?

For nodes which are not connected to regular shell elements, and therefore do not have director information assigned, these methods all return a zero-valued vector.

If not otherwise specified, the current and rest director values are computed automatically from the surrounding (regular) shell elements, with their value {\bf d} being computed from

{\bf d}=\sum t_{i}\,{\bf n}_{i}

where t_{i} is the value of the defaultThickness property and {\bf n}_{i} is the surface normal of the i-th surrounding regular shell element. However, if necessary, it is also possible to explicitly assign these values, using the methods

  setDirector (Vector3d dir);     // set the current director
  setRestDirector (Vector3d dir); // set the rest director

ArtiSynth FEM nodes can currently support only one director, which is shared by all regular shell elements associated with that node. This effectively means that all such elements must belong to the same “surface”, and that two intersecting surfaces cannot share the same nodes.

As indicated above, shell elements can also be instantiated as membrane elements, which do not exhibit bending stiffness and therefore do not require director information. The regular/membrane distinction is specified in the element’s constructor. For example, ShellTriElement and ShellQuadElement each have constructors with the signatures:

  ShellTriElement (FemNode3d n0, FemNode3d n1, FemNode3d n2,
    double thickness, boolean membrane);
  ShellQuadElement (FemNode3d n0, FemNode3d n1, FemNode3d n2, FemNode3d n3,
    double thickness, boolean membrane);

The thickness argument specifies the defaultThickness property, while membrane determines whether or not the element is a membrane element.

While membrane elements do not require explicit director information stored at the nodes, they do make use of an inferred director that is parallel to the element’s surface normal, and has a constant length equal to the element’s defaultThickness property. This gives the element a virtual volume, which (as with regular elements) is used to determine 3D strains and to compute the element’s mass from it’s density.

6.1.2.4 Meshes


The geometry associated with a finite element model consists of a collection of meshes (e.g. PolygonalMesh, PolylineMesh, PointMesh) that move along with the model in a way that maintains the shape function interpolation equation (6.2) at each vertex location. These geometries can be used for visualizations, or for physical interactions like collisions. However, they have no physical properties themselves. FEM geometries will be discussed in more detail in Section 6.3. The list of meshes can be obtained with the method

  MeshComponentList<FemMeshComp> getMeshComps();

6.1.3 Materials

The stress-strain relationship within each element is defined by a “material” delegate object, implemented by a subclass of FemMaterial. This material object is responsible for implementing the functions

   void computeStressAndTangent (...)

which computes the stress tensor and (optionally) the tangent stiffness matrix at each integration point, based on the current local deformation at that point.

The default material type is LinearMaterial, where stress is related to strain through:

\displaystyle\sigma({\bf x})=D\,\epsilon({\bf x}), (6.3)
\displaystyle\text{where }\quad D=\begin{bmatrix}\lambda+2\mu&\lambda&\lambda&%
0&0&0\\
\lambda&\lambda+2\mu&\lambda&0&0&0\\
\lambda&\lambda&\lambda+2\mu&0&0&0\\
0&0&0&\mu&0&0\\
0&0&0&0&\mu&0\\
0&0&0&0&0&\mu\end{bmatrix},\quad\lambda=\dfrac{E\nu}{(1+\nu)(1-2\nu)},\quad\mu%
=\dfrac{E}{2(1+\nu)},

\sigma is the standard 6\times 1 stress vector, \epsilon is the strain vector, E is the Young’s Modulus, and \nu is Poisson’s ratio. This linear material uses a corotational formulation, so rotations are removed per element before computing the strain [13]. To enable or disable this corotational formulation, use LinearMaterial.setCorotated(boolean).

All material models, including linear and nonlinear, are available in the package artisynth.core.materials. A list of common materials is provided in Table 6.3. Those that are subclasses of IncompressibleMaterial allow for incompressibility.

Table 6.3: Commonly used FEM materials
Material Parameters
LinearMaterial E Young’s modulus
\nu Poisson’s ratio
corotated corotational formulation
StVenantKirchoffMaterial E Young’s modulus
\nu Poisson’s ratio
NeoHookeanMaterial E Young’s modulus
\nu Poisson’s ratio
IncompNeoHookeanMaterial G shear modulus
\kappa bulk modulus
MooneyRivlinMaterial C_{10},C_{01},C_{20},C_{02} distortional parameters
\kappa bulk modulus
OgdenMaterial \mu_{1},\ldots,\mu_{6} material parameters
\alpha_{1},\ldots,\alpha_{6}
\kappa bulk modulus

6.1.4 Boundary conditions

Boundary conditions can be implemented in one of several ways:

  1. 1.

    Explicitly setting FEM node positions/velocities

  2. 2.

    Attaching FEM nodes to other dynamic components

  3. 3.

    Enabling collisions

To enforce an explicit (Dirichlet) boundary condition for a set of nodes, their dynamic property must be set to false. This notifies ArtiSynth that the state of these nodes (both position and velocity) will be controlled parametrically. By disabling dynamics, a fixed boundary condition is applied. For a moving boundary, positions and velocities of the boundary nodes must be explicitly set every timestep. This can be accomplished with either a Controller (see Section 5.3) or an InputProbe (see Section 5.4). Note that both the position and velocity of the nodes should be explicitly set for consistency.

Another type of supported boundary condition is to attach FEM nodes to other components, including particles, springs, rigid bodies, and locations within other FEM elements. Here, the node is still considered dynamic, but its motion is coupled to that of the attached component through a constraint mechanism. Attachments will be discussed further in Section 6.4.

Finally, the boundary of an FEM can be constrained by enabling collisions with other components. This will be covered in Section 6.11.