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 [5].
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
void setDensity (double density) |
Sets the density. |
double getDensity() |
Gets the density. |
void setMaterial (FemMaterial mat) |
Sets the FEM’s material. |
FemMaterial getMaterial() |
Gets the FEM’s material. |
void setParticleDamping (double d) |
Sets the particle (mass) damping. |
double getParticleDamping() |
Gets the particle (mass) damping. |
void setStiffnessDamping (double d) |
Sets the stiffness damping. |
double getStiffnessDamping() |
Gets the stiffness damping. |
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 individual 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
(6.0) |
where is the value of particleDamping, is the value of stiffnessDamping, is the FEM model’s lumped mass matrix, is the FEM’s stiffness matrix, and 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”.
Each FemModel3d contains several lists of subcomponents:
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).
The volumetric model elements. These define the 3D sub-units over which the system is numerically integrated.
The shell elements. These define additional 2D sub-units over which the system is numerically integrated.
The geometry in the model. This includes the surface mesh, and any other embedded geometries.
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.
Optional field components which can be used to interpolate application-defined quantities over the FEM model’s domain. Fields are described in detail in Chapter 7.
The nodes, elements and meshes components are illustrated in Figure 6.2.
(a) FEM model | (b) Nodes | (c) Elements | (d) Geometry |
The set of nodes belong to a finite element model can be obtained by the method
PointList<FemNode3d> getNodes() |
Returns the 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).
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:
(6.1) |
where is the displacement of the th node that is adjacent to the element, and 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 |
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() |
Returns the list of volumetric elements. |
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.
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:
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() |
Returns the list of shell elements. |
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() |
Returns a list of all elements. |
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 [11], and more specifically the front/back node formulation [12]. 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.
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
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 being computed from
where is the value of the defaultThickness property and is the surface normal of the -th surrounding regular shell element. However, if necessary, it is also possible to explicitly assign these values, using the methods
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:
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.
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.1) 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
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
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:
(6.2) | |||
is the standard stress vector, is the strain vector, is the Young’s Modulus, and is Poisson’s ratio. This linear material uses a corotational formulation, so rotations are removed per element before computing the strain [16]. 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.
Material | Parameters | |
---|---|---|
LinearMaterial | Young’s modulus | |
Poisson’s ratio | ||
corotated | corotational formulation | |
StVenantKirchoffMaterial | Young’s modulus | |
Poisson’s ratio | ||
NeoHookeanMaterial | Young’s modulus | |
Poisson’s ratio | ||
IncompNeoHookeanMaterial | shear modulus | |
bulk modulus | ||
MooneyRivlinMaterial | distortional parameters | |
bulk modulus | ||
OgdenMaterial | material parameters | |
bulk modulus |
Boundary conditions can be implemented in one of several ways:
Explicitly setting FEM node positions/velocities
Attaching FEM nodes to other dynamic components
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 Chapter 8.