Finite Element Models

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:

These properties can be set and retrieved using the methods

Keep in mind that ArtiSynth is essentially “unitless” (Section LABEL:sec:mechii:units), 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 Sections 1.1.3 and 1.10. 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 1.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

(1.1)

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 (LABEL:eqn:pointdamping), which is why it is referred to as “particle damping”.

Each FemModel3d contains several lists of subcomponents:

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 1.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 Chapter LABEL:sec:fields.

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

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.

All supported materials are presented in detail in Section 1.10. The default material type is LinearMaterial, which linearly maps Cauchy strain onto Cauchy stress via

where is the elasticity tensor determined from Young’s modulus and Poisson’s ratio. Anisotropic linear materials are provided by TransverseLinearMaterial and AnisotropicLinearMaterial. Linear materials in ArtiSynth implement corotation, which removes the rotation from the deformation gradient and allows Cauchy strain to be applied to large deformations [ngan:fem:2008, muller2004interactive]. Nonlinear hyperelastic materials are also supported, including MooneyRivlinMaterial, OgdenMaterial, YeohMaterial, ArrudaBoyceMaterial, and VerondaWestmannMaterial.

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 LABEL:ControllersAndMonitors:sec) or an InputProbe (see Section LABEL:Probes:sec). 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 1.4.

Finally, the boundary of an FEM can be constrained by enabling collisions with other components. This will be covered in Chapter LABEL:ContactAndCollision:sec.

For simple shapes such as beams and ellipsoids, there are factory methods to automatically build the node and element structure. These methods are found in the FemFactory class. Some common methods are

The inputs specify the dimensions, resolution, and potentially the type of element to use. The following code creates a basic beam made up of hexahedral elements:

For more complex geometries, volumetric meshes can be loaded from external files. A list of supported file types is provided in Table 1.3. To load a geometry, an appropriate file reader must be created. Readers capable of reading FEM models implement the interface FemReader, which has the method

Additionally, many FemReader classes have static methods to handle the loading of files for convenience.

The following code snippet demonstrates how to load a model using the AnsysReader.

The method PathFinder.getSourceRelativePath() is used to find a path within the ArtiSynth source tree that is relative to the current model’s source file (Section LABEL:PathFinder:sec). Note the try-catch block. Most of these readers throw an IOException if the read fails.

There are two ways an FEM model can be generated from a surface: by using a FEM mesh generator, and by extruding a surface along its normal direction.

ArtiSynth has the ability to interface directly with the TetGen library (http://tetgen.org) to create a tetrahedral volumetric mesh given a closed and manifold surface. The main Java class for calling TetGen directly is TetgenTessellator. The tessellator has several advanced options, allowing for the computation of convex hulls, and for adding points to a volumetric mesh. For simply creating an FEM from a surface, there is a convenience routine within FemFactory that handles both mesh generation and constructing a FemModel3d:

If quality , then points will be added in an attempt to bound the maximum radius-edge ratio (see the -q switch for TetGen). According to the TetGen documentation, the algorithm usually succeeds for a quality ratio of 1.2.

It’s also possible to create thin layer of elements by extruding a surface along its normal direction.

For example, to create a two-layer slice of elements centered about a surface of a tendon mesh, one might use

For this type of extrusion, triangular faces become wedge elements, and quadrilateral faces become hexahedral elements.

Note: for extrusions, no care is taken to ensure element quality; if the surface has a high curvature relative to the total extrusion thickness, then some elements will be inverted.

A finite element model’s structure can also be manually constructed in code. FemModel3d has the methods:

For an element to successfully be added, all its nodes must already have been added to the model. Nodes can be constructed from a 3D location, and elements from an array of nodes. A convenience routine is available in FemElement3d that automatically creates the appropriate element type given the number of nodes (Table 1.1):

Be aware of node orderings when supplying nodes. For linear elements, ArtiSynth uses a clockwise convention with respect to the outward normal for the first face, followed by the opposite node(s). To determine the correct ordering for a particular element, check the coordinates returned by the function FemElement3dBase.getNodeCoords(). This returns the concatenated coordinate list for an “ideal” element of the given type.

A complete application model that implements a simple FEM beam is given below.

This example can be found in artisynth.demos.tutorial.FemBeam. The build() method first creates a MechModel and FemModel3d. A FEM beam is created using a factory method on line 36. This beam is centered at the origin, so its length extends from -length/2 to length/2. The density, damping and material properties are then assigned.

On lines 45–49, a fixed boundary condition is set to the left-hand side of the beam by setting the corresponding nodes to be non-dynamic. Due to numerical precision, a small EPSILON buffer is required to ensure all left-hand boundary nodes are identified (line 46).

Rendering properties are then assigned to the FEM model on line 52. These will be discussed further in Section 1.12.

By default, every FemModel3d has an auto-generated geometry representing the “surface mesh”. The surface mesh consists of all un-shared element faces (i.e. the faces of individual elements that are exposed to the world), and its vertices correspond to the nodes that make up those faces. As the FEM nodes move, so do the mesh vertices due to the attachment framework.

The surface mesh can be obtained using one of the following functions in FemModel3d:

The first returns the surface complete with attachment information. The latter method directly returns the PolygonalMesh that is controlled by the FEM.

It is possible to manually set the surface mesh:

However, doing so is normally not necessary. It is always possible to add additional mesh geometries to a finite element model, and the visibility settings can be changed so that the default surface mesh is not rendered.

Any geometry of type MeshBase can be added to a FemModel3d. To do so, first position the mesh so that its vertices are in the desired locations inside the FEM, then call one of the FemModel3d methods:

The latter is a convenience routine that also gives the newly embedded FemMeshComp a name.

A complete model demonstrating embedding a mesh is given below.

This example can be found in artisynth.demos.tutorial.FemEmbeddedSphere. The model is very similar to FemBeam. A MechModel and FemModel3d are created and added. At line 41, a PolygonalMesh of a sphere is created using a factory method. The sphere is already centered inside the beam, so it does not need to be repositioned. At Line 42, the sphere is embedded inside model fem, creating a FemMeshComp with the name “sphere”. The full model is shown in Figure 1.5.

Since FemNode3d is a subclass of Particle, the same methods described in Section LABEL:sec:mech:pointattachments for attaching particles to other particles and frames are available. For example, we can attach an FEM node to a rigid body using a either a statement of the form

or the following equivalent statement which does the same thing:

Both of these create a PointFrameAttachment between a rigid body (called body) and an FEM node (called node) and then adds it to the MechModel mech.

One can also attach the nodes of one FEM model to the nodes of another using statements like

or

which attaches node2 to node1.

The following model demonstrates attaching an FEM beam to a rigid block.

This model extends the FemBeam example of Section 1.2.5. The build() method then creates and adds a RigidBody block (lines 18–20). On line 21, the block is repositioned to the side of the beam to prepare for the attachment. On lines 24–28, all right-most nodes of the beam are then set to be attached to the block using a PointFrameAttachment. In this case, the attachments are explicitly created. They could also have been attached using

Typically, nodes do not align in a way that makes it possible to connect them to other FEM models and/or points based on simple point-to-node attachments. Instead, we use a different mechanism that allows us to attach a point to an arbitrary location within an FEM model. This is done using an attachment component of type PointFem3dAttachment, which implements an attachment where the position and velocity of the attached point is determined by a weighted sum of the positions and velocities of selected fem nodes:

(1.3)

Any force acting on the attached point is then propagated back to the nodes, according to the relation

(1.4)

where is the force acting on node due to . This relation can be derived based on the conservation of energy. If is embedded within a single element, then the are simply the element nodes and the are corresponding shape function values; this is known as an element-based attachment. On the other hand, as described below, it is sometimes desirable to form an attachment using a more general set of nodes that extends beyond a single element; this is known as a nodal-based attachment (Section 1.4.8).

An element-based attachment can be created using a code fragment of the form

First, a PointFem3dAttachment is created for the point pnt. Next, setFromElement() is used to determine the nodal weights within the element elem for the specified position (which in this case is simply the point’s current position). To do this, it computes the “natural coordinates” coordinates of the position within the element. For this to be guaranteed to work, the position should be on or inside the element. If natural coordinates cannot be found, the method will return false and the nearest estimates coordinates will be used instead. However, it is sometimes possible to find natural coordinates outside a given element as long as the shape functions are well-defined. Finally, the attachment is added to the model.

More conveniently, the exact same functionality is provided by the attachPoint() method in MechModel:

This creates an attachment identical to that created by the previous code fragment.

Often, one does not want to have to determine the element to which a point should be attached. In that case, one can call

or, equivalently,

This will find the nearest element to the node in question and use that to create the attachment. If the node is outside the FEM model, then it will be attached to the nearest point on the FEM’s surface.

The following model demonstrates how to attach two FEM models together:

This example can be found in artisynth.demos.tutorial.FemBeamWithFemSphere. The model extends FemBeam, adding a finite element sphere and coupling them together. The sphere is created and added on lines 18–28. It is read from TetGen-generated files using the TetGenReader class. The model is then scaled to match the dimensions of the current model, and transformed to the right side of the beam. To create attachments, the code first checks for any nodes that belong to the sphere that fall inside the beam using the FemModel3d.findContainingElement(Point3d) method (line 36), which returns the containing element if the point is inside the model, or null if the point is outside. Internally, this spatial search uses a bounding volume hierarchy for efficiency (see BVTree and BVFeatureQuery). If the point is contained within the beam, then mech.attachPoint() is used to attach it to the nodes of the element (line 39).

While it is straightforward to connect nodes to rigid bodies or other FEM nodes or elements, it is often necessary to determine which nodes to attach. This was evident in the example of Section 1.4.4, which attached nodes of an FEM sphere that were found to be inside an FEM beam.

As in that example, finding the nodes to attach can often be done using geometric queries. For example, we often select nodes based on how close they are to the other body we wish to attach to.

Various prolixity queries are available for this task. To find the distance of a point to a polygonal mesh, we can use the following PolygonalMesh methods,

where the latter method returns 1 if the point is inside and 0 otherwise. For checking the distance of an FEM node, pnt can be obtained from node.getPosition() (or possibly node.getRestPosition()). For example, to find all nodes within a distance tol of the surface of a rigid body, we could use the code fragment:

If we want to check only nodes that lie on the FEM surface, then we can filter them using the FemModel3d method isSurfaceNode():

Most of the mesh-based query methods work only for triangular meshes. The PolygonalMesh method isTriangular() can be used to determine if the mesh is triangular. If it is not, it can made triangular by calling triangulate(), although in general this should be done during model construction before the mesh-based component has been added to the model.

For connecting an FEM model to another FEM model, FemModel3d provides a number of query methods:

All the above queries are based on the FEM model’s current nodal positions. The method
findNearestElement(near,pnt,filter) allows the application to specify a FemModel.ElementFilter to restrict the elements that are searched.

The argument near that appears in some of the queries is an optional argument which, if not null, returns the location of the corresponding nearest point on the element. The distance from pnt to the element can then be found using

   near.distance (pnt);

If the resulting distance is 0, then the point is on or inside the element. Otherwise, the point is outside the element, and if no element filters were used in the query, outside the FEM model itself.

Typically, it is preferred attach a point to an element only if it lies on or inside an element. However, it is possible to attach points outside an element as long as the system is able to determine appropriate element “coordinates” for that point (which it may not be able to do if the point is far away). In addition, the motion behavior of an exterior attached point may sometimes appear counterintuitive.

The FemModel3d element and node queries can be used in a variety of ways.

findNearestNodes() can be used to find all nodes within a certain distance of a point, as part of the process of making nodal-based attachments (Section 1.4.8).

findNearestNode() is used in the FemMuscleBeam example (Section 1.9.5) to determine if a desired muscle point is near enough to a node to use that node directly, or if a marker should be created.

As another example, suppose we wish to connect the surface nodes of an FEM model femA to the surface elements of another model femB if they lie within a prescribed distance tol of the surface of femB. Then we could use the following code:

Finally, it is possible to identify nodes on the surface of an FEM model according to whether they belong to specific features, such as a smooth patch or a sharp edge line. Methods for doing this are provided as static methods in the class FemQuery, and include:

Details of how these methods work are given in their API documentation. They use the notion of a bend angle, which is the absolute value of the angle between two faces about their common edge. A patch is defined by a collection of faces whose bend angles do not exceed a minimum value, while an edge line is a collection of edges with bend angles not below a maximum value. The feature methods start with an initial node (node0) and then grow the requested feature out from there. For example, suppose we have a regular hexahedral FEM grid, and we wish to find all the nodes on one of the faces. If we know one of the nodes on the face, then we can find all of the nodes using findPatchNodes:

Note that the feature query above uses a maximum bend angle of . Because grid faces are flat, this choice is somewhat arbitrary; any angle larger than 0 (within machine precision) would also work.

Often, it is most convenient to select nodes in the ArtiSynth viewer. For this, a node selection tool is available (Figure 1.8), as described in the section “Selecting FEM nodes” of the ArtiSynth User Interface Guide). It allows nodes to be selected in various ways, including the usual click, drag box and elliptical selection tools, as well as specialized operations that select nodes based on patches, edge lines, and distances to other bodies. Once nodes have been selected, their numbers can be saved to a node number file which can then be read by a model’s build() method to determine which nodes to connect to some body.

Node number files are text files with a very simple format, consisting of integers (the node numbers), separated by white space. There is no limit to how many numbers can appear on a line but typically this is limited to ten or so to make the file more readable. Optionally, the numbers can be surrounded by square brackets ([ ]). The special character ‘#’ is a comment character, commenting out all characters from itself to the end of the current line. For a file containing the node numbers 2, 12, 4, 8, 23 and 47, the following formats are all valid:

Once node numbers have been identified in the viewer and saved in a file, they can be read by the build() method using a NodeNumberReader. For convenience, this class supplies two static methods for extracting the FEM nodes specified by the numbers in a file:

Extracted nodes can be used to set boundary conditions or form connections with other bodies. For example, suppose we wish to connect a face of an FEM model fem to a rigid body block, using a set of nodes specified in a file called faceNodes.txt. A code fragment to accomplish this could be the following:

The process of selecting nodes in the viewer, saving them in a file, and using them in a model can be done iteratively: if the selected nodes need to be adjusted, one can reopen the node selection tool, load the selection from file, adjust the selection, and resave the file, without needing to make any modifications to the model’s build() code.

If desired, one can also determine a set of nodes in code, and then write their numbers to a file using the class NodeNumberWriter, which supplies static methods for writing number files:

For example:

The LumbarFrameSpring example in Section LABEL:LumbarFrameSpring:sec uses a frame spring to connect two simplified lumbar vertebrae. However, it is also possible to use an FEM model in place of a frame spring, possibly providing a more realistic model of the intervertebral disc. A simple model which does this is defined in

  artisynth.demos.tutorial.LumbarFEMDisk

The initial source code is similar to that for LumbarFrameSpring, but differs in the section where the FEM disk replaces the FrameSpring:

The simplified FEM model representing the “disk” is created at lines 57-61, using a torus-shaped model created by FemFactory. It is then repositioning using transformGeometry() ( Section LABEL:TransformingGeometry:sec) to place it between the vertebrae (line 64-66). After the FEM model is positioned, we find which nodes are within a distance tol of each vertebral surface and attach them to the appropriate body (lines 69-81).

To run this example in ArtiSynth, select All demos > tutorial > LumbarFEMDisk from the Models menu. The model should load and initially appear as in Figure 1.9. The behavior is best seem by running the model and using the pull controller to exert forces on the upper vertebra.

The example of Section 1.4.4 uses element-based attachments to connect the nodes of one FEM to elements of another. As mentioned above, element-based attachments assume that the attached point is associated with a specific FEM model element. While this often gives good results, there are situations where it may be desirable to distribute the connection more broadly among a larger set of nodes.

In particular, this is sometimes the case when connecting FEM models to point-to-point springs. The end-point of such a spring may end up exerting a large force on the FEM, and then if the number of nodes to which the end-point is attached are too small, the resulting forces on these nodes (Equation 1.4) may end up being too large. In other words, it may be desirable to distribute the spring’s force more evenly throughout the FEM model.

To handle such situations, it is possible to create a nodal-based attachment in which the nodes and weights are explicitly specified. This involves explicitly creating a PointFem3dAttachment for the point or particle to be attached and the specifying the nodes and weights directly,

where nodes and weights are arrays of FemNode and double, respectively. It is up to the application to determine these.

PointFem3dAttachment provides several methods for explicitly specifying nodes and weights. The signatures for these include:

The last two methods determine the weights automatically, using an inverse-distance-based calculation in which each weight is initially computed as

(1.4)

where is the distance from node to pos and is the maximum distance. The weights are then adjusted to ensure that they sum to one and that the weighted sum of the nodes equals pos. In some cases, the specified nodes may not provide enough support for the last condition to be met, in which case the methods return false.

The model demonstrating the difference between element and nodal-based attachments is defined in

  artisynth.demos.tutorial.PointFemAttachment

It creates two FEM models, each with a single point-to-point spring attached to a particle at their center. The model at the top (fem1 in the code below) is connected to the particle using an element-based attachment, while the lower model (fem2 in the code) is connected using a nodal-based attachment with a larger number of nodes. Figure 1.10 shows the result after the model is run until stable. The element-based attachment results in significantly higher deformation in the immediate vicinity around the attachment, while for the nodal-based attachment, the deformation is distributed much more evenly through the model.

The build method and some of the auxiliary code for this model is shown below. Code for the other auxiliary methods, including addFem(), addParticle(), addSpring(), and setAttachedNodesWhite(), can be found in the actual source file.

The build() method begins by creating a MechModel and then adding to it two FEM beams (created using the auxiliary method addFem(). Rendering of each FEM model’s surface is then set up to show strain values (setSurfaceRendering(), lines 41 and 43). The surface meshes themselves are also redefined to exclude the frontmost elements, allowing the strain values to be displayed closer model centers. This redefinition is done using calls to createSurfaceMesh() (lines 40, 41) with a custom ElementFilter defined at lines 3-12.

Next, the end-point particles for the axial springs are created (using the auxiliary method addParticle(), lines 46-49), and particle m1 is attached to fem1 using mech.attachPoint() (line 52), which creates an element-based attachment at the point’s current location. Point m2 is then attached to fem2 using a nodal-based attachment. The nodes for these are collected as the union of all nodes for a specified set of elements (lines 58-59, and the method collectNodes() defined at lines 16-25). These are then used to create a nodal-based attachment (lines 61-63), where the weights are determined automatically using the method associated with equation (1.4).

Finally, the springs are created (auxiliary method addSpring(), lines 66-67), the nodes associated for each attachment are set to render as white spheres (setAttachedNodesWhites(), lines 70-71), and the particles are set to render as green spheres.

To run this example in ArtiSynth, select All demos > tutorial > PointFemAttachment from the Models menu. Running the model will cause it to settle into the state shown in Figure 1.10. Selecting and dragging one of the spring anchor points at the right will cause the spring tension to vary and further illustrate the difference between the element and nodal-based attachments.

A complete application model that employs a fem marker as an anchor for a spring is given below.

This example can be found in artisynth.demos.tutorial.FemBeamWithMuscle. This model extends the FemBeam example, adding a FemMarker for the spring to attach to. The method createMarker(...) on lines 29–35 is used to create and add a marker to the FEM. Since the element is initially set to null, when it is added to the FEM, the model searches for the containing or nearest element. The loaded model is shown in Figure 1.11.

A model illustrating how to connect frames to an FEM model is defined in

  artisynth.demos.tutorial.FrameFemAttachment

It creates an FEM beam, along with a rigid body block and a massless coordinate frame, that are then attached to the beam using nodal and element-based attachments. The build method is shown below:

Lines 3-22 create a MechModel and populate it with an FEM beam and a rigid body box. Next, a basic Frame is created, with a specified pose and an axis length of 0.3 (to allow it to be seen), and added to the MechModel (lines 25-28). It is then attached to the FEM beam using an element-based attachment (line 30). Meanwhile, the box is attached to using a nodal-based attachment, created from all the nodes associated with elements 22 and 31 (lines 33-36). Finally, all attachment nodes are set to be rendered as green spheres (lines 39-41).

To run this example in ArtiSynth, select All demos > tutorial > FrameFemAttachment from the Models menu. Running the model will cause it to settle into the state shown in Figure 1.12. Forces can interactively be applied to the attached block and frame using the pull tool, causing the FEM model to deform (see the section “Pull Manipulation” in the ArtiSynth User Interface Guide).

The ability to connect frames to FEM models, as described in Section 1.6, makes it possible to interconnect different FEM models directly using joints, as described in Section LABEL:JointsAndConnectors:sec. This is done internally by using FrameFem3dAttachments to connect frames C and D of the joint (Figure LABEL:jointExample:fig) to their respective FEM models.

As indicated in Section LABEL:CreatingJoints:sec, most joints have a constructor of the form

that creates a joint connecting bodyA to bodyB, with the initial pose of the D frame given (in world coordinates) by TDW. The same body and transform settings can be made on an existing joint using the method setBodies(bodyA, bodyB, TDW). For these constructors and methods, it is possible to specify FEM models for bodyA and/or bodyB. Internally, the joint then creates a FrameFem3dAttachment to connect frame C and/or D of the joint (See Figure LABEL:jointExample:fig) to the corresponding FEM model.

However, unlike joints involving rigid bodies or frames, there are no associated or transforms (since there is no fixed frame within an FEM to define such transforms). Methods or constructors which utilize or can therefore not be used with FEM models.

A model connecting two FEM beams by a joint is defined in

  artisynth.demos.tutorial.JointedFemBeams

It creates two FEM beams and connects them via a special slotted-revolute joint. The build method is shown below: