6 Finite Element Models

6.4 Connecting FEM models to other components

To couple FEM models to other dynamic components, the “attachment” mechanism described in Section 1.2 is used. This involves creating and adding to the model attachment components, which are instances of DynamicAttachment, as described in Section 3.6. Common point-based attachment classes are listed in Table 6.5.

Table 6.5: Point-based attachments
Attachment Description
PointParticleAttachment Attaches one “point” to one “particle”
PointFrameAttachment Attaches one “point” to one “frame”
PointFem3dAttachment Attaches one “point” to a linear combination of FEM nodes

FEM models are connected to other model components by attaching their nodes to various components. This can be done by creating an attachment object of the appropriate type, and then adding it to the MechModel using

  addAttachment (DynamicAttachment attach); // adds an attachment constraint

There are also convenience routines inside MechModel that will create the appropriate attachments automatically (see Section 3.6.1).

All attachments described in this section are based around FEM nodes. However, it is also possible to attach frame-based components (such as rigid bodies) directly to an FEM, as described in Section 6.6.

6.4.1 Connecting nodes to rigid bodies or particles

Since FemNode3d is a subclass of Particle, the same methods described in Section 3.6.1 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

   mech.addAttachment (new PointFrameAttachment(body, node));

or the following equivalent statement which does the same thing:

   mech.attachPoint (node, body);

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

   mech.addAttachment (new PointParticle (node1, node2));

or

   mech.attachPoint (node2, node1);

which attaches node2 to node1.

6.4.2 Example: connecting a beam to a block

Figure 6.6: FemBeamWithBlock model loaded into artisynth.

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

1 package artisynth.demos.tutorial;
2
3 import java.io.IOException;
4
5 import maspack.matrix.RigidTransform3d;
6 import artisynth.core.femmodels.FemNode3d;
7 import artisynth.core.mechmodels.PointFrameAttachment;
8 import artisynth.core.mechmodels.RigidBody;
9
10 public class FemBeamWithBlock extends FemBeam {
11
12    public void build (String[] args) throws IOException {
13
14       // Build simple FemBeam
15       super.build (args);
16
17       // Create a rigid block and move to the side of FEM
18       RigidBody block = RigidBody.createBox (
19          "block", width/2, 1.2*width, 1.2*width, 2*density);
20       mech.addRigidBody (block);
21       block.setPose (new RigidTransform3d (length/2+width/4, 0, 0));
22
23       // Attach right-side nodes to rigid block
24       for (FemNode3d node : fem.getNodes()) {
25          if (node.getPosition().x >= length/2-EPS) {
26             mech.addAttachment (new PointFrameAttachment (block, node));
27          }
28       }
29    }
30
31 }

This model extends the FemBeam example of Section 6.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

   mech.attachPoint (node, block);  // attach node to rigid block

6.4.3 Connecting nodes directly to elements

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 {\bf p} and velocity {\bf u} of the attached point is determined by a weighted sum of the positions {\bf x}_{k} and velocities {\bf u}_{k} of selected fem nodes:

{\bf p}=\sum\alpha_{k}{\bf x}_{k} (6.4)

Any force {\bf f} acting on the attached point is then propagated back to the nodes, according to the relation

{\bf f}_{k}=\alpha_{k}{\bf f} (6.5)

where {\bf f}_{k} is the force acting on node k due to {\bf f}. This relation can be derived based on the conservation of energy. If {\bf p} is embedded within a single element, then the {\bf x}_{k} are simply the element nodes and the \alpha_{i} 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 6.4.7).

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

   PointFem3dAttachment ax = new PointFem3dAttachment(pnt);
   ax.setFromElement (pnt.getPosition(), elem);
   mech.addAttachment (ax);

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:

   mech.attachPoint (pnt, elem);

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

   PointFem3dAttachment ax = new PointFem3dAttachment(pnt);
   ax.setFromFem (pnt.getPosition(), fem);
   mech.addAttachment (ax);

or, equivalently,

   mech.attachPoint (pnt, fem);

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.

6.4.4 Example: connecting two FEMs together

Figure 6.7: FemBeamWithFemSphere model loaded into ArtiSynth.

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

1 package artisynth.demos.tutorial;
2
3 import java.io.IOException;
4
5 import maspack.matrix.RigidTransform3d;
6 import artisynth.core.femmodels.*;
7 import artisynth.core.materials.LinearMaterial;
8 import maspack.util.PathFinder;
9
10 public class FemBeamWithFemSphere extends FemBeam {
11
12    public void build (String[] args) throws IOException {
13
14       // Build simple FemBeam
15       super.build (args);
16
17       // Create a FEM sphere
18       FemModel3d femSphere = new FemModel3d("sphere");
19       mech.addModel(femSphere);
20       // Read from TetGen file
21       TetGenReader.read(femSphere,
22          PathFinder.getSourceRelativePath(FemModel3d.class, "meshes/sphere2.1.node"),
23          PathFinder.getSourceRelativePath(FemModel3d.class, "meshes/sphere2.1.ele"));
24       femSphere.scaleDistance(0.22);
25       // FEM properties
26       femSphere.setDensity(10);
27       femSphere.setParticleDamping(2);
28       femSphere.setMaterial(new LinearMaterial(4000, 0.33));
29
30       // Reposition FEM to side of beam
31       femSphere.transformGeometry( new RigidTransform3d(length/2+width/2, 0, 0) );
32
33       // Attach sphere nodes that are inside beam
34       for (FemNode3d node : femSphere.getNodes()) {
35          // Find element containing node (if exists)
36          FemElement3d elem = fem.findContainingElement(node.getPosition());
37          // Add attachment if node is inside "fem"
38          if (elem != null) {
39             mech.attachPoint(node, elem);
40          }
41       }
42
43       // Set render properties
44       setRenderProps(femSphere);
45
46    }
47
48 }

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).

6.4.5 Finding which nodes to attach

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 6.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 promixity queries are available for this task. To find the distance of a point to a polygonal mesh, we can use the following PolygonalMesh methods,

   double distanceToPoint (Point3d pnt)
   int isInside (Point3d pnt) // for closed meshes only

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:

   RigidBody body;
   FemModel3d fem;
   ...
   double tol = 0.001;
   PolygonalMesh surface = body.getSurfaceMesh();
   ArrayList<FemNode3d> nearNodes = new ArrayList<FemNode3d>();
   for (FemNode3d n : fem.getNodes()) {
      if (surface.distanceToPoint (n.getPosition()) < tol) {
         nearNodes.add (n);
      }
   }

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

   for (FemNode3d n : fem.getNodes()) {
      if (fem.isSurfaceNode (n) &&
          surface.distanceToPoint (n.getPosition()) < tol) {
         nearNodes.add (n);
      }
   }

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:

   /* nearest element queries */
   // find the nearest element (shell or volumetric) to pnt:
   FemElement3dBase findNearestElement (Point3d near, Point3d pnt)
   FemElement3dBase findNearestElement (
      Point3d near, Point3d pnt, ElementFilter filter)
   // find the nearest surface element (shell or volumetric) to pnt:
   FemElement3dBase findNearestSurfaceElement(Point3d near, Point3d pnt)
   // find the nearest volumetric element to pnt:
   FemElement3d findNearestVolumetricElement (Point3d near, Point3d pnt)
   // find the nearest shell element to pnt:
   ShellElement3d findNearestShellElement (Point3d near, Point3d pnt)
   // find the volumetric element (if any) containing pnt:
   FemElement3d findContainingElement(Point3d pnt)
   /* nearest node queries */
   // find the nearest node to pnt that is within maxDist:
   FemNode3d findNearestNode (Point3d pnt, double maxDist)
   // find all nodes that are within maxDist of pnt:
   ArrayList<FemNode3d> findNearestNodes (Point3d pnt, double maxDist)

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 6.4.7).

findNearestNode() is used in the FemMuscleBeam example (Section 6.9.4) 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:

   MechModel mech;
   FemModel3d femA;
   FemModel3d femB;
   ...
   double tol = 0.001;
   Point3d near = new Point3d();
   for (FemNode3d n : femA.getNodes()) {
      if (femA.isSurfaceNode(n)) {
         FemElement3dBase elem =
            femB.findNearestSurfaceElement (near, n.getPosition());
         if (elem != null && near.distance(n.getPosition()) <= tol) {
            // attach if within distance
            mech.attachPoint (n, elem);
         }
      }
   }
Figure 6.8: LumbarFEMDisk loaded into ArtiSynth, showing a simplified FEM model connecting two vertebrae.

6.4.6 Example: two bodies connected by an FEM “spring”

The LumbarFrameSpring example in Section 3.5.4 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:

56       // create a torus shaped FEM model for the disk
57       FemModel3d fem = new FemModel3d();
58       fem.setDensity (1500);
59       fem.setMaterial (new LinearMaterial (20000, 0.4));
60       FemFactory.createHexTorus (fem, 0.011, 0.003, 0.008, 16, 30, 2);
61       mech.addModel (fem);
62
63       // position it betweem the disks
64       double DTOR = Math.PI/180.0;
65       fem.transformGeometry (
66          new RigidTransform3d (-0.012, 0.0, 0.040, 0, -DTOR*25, DTOR*90));
67
68       // find and attach nearest nodes to either the top or bottom mesh
69       double tol = 0.001;
70       for (FemNode3d n : fem.getNodes()) {
71          // top vertebra
72          double d = lumbar1.getSurfaceMesh().distanceToPoint(n.getPosition());
73          if (d < tol) {
74             mech.attachPoint (n, lumbar1);
75          }
76          // bottom vertebra
77          d = lumbar2.getSurfaceMesh().distanceToPoint(n.getPosition());
78          if (d < tol) {
79             mech.attachPoint (n, lumbar2);
80          }
81       }
82       // set render properties ...
83       fem.setSurfaceRendering (SurfaceRender.Shaded);
84       RenderProps.setFaceColor (fem, new Color (153/255f, 153/255f, 1f));
85       RenderProps.setFaceColor (mech, new Color (238, 232, 170)); // bone color
86    }
87 }

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 4.8) 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 6.8. The behavior is best seem by running the model and using the pull controller to exert forces on the upper vertebra.

6.4.7 Nodal-based attachments

The example of Section 6.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 6.5) 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,

   PointFem3dAttachment ax = new PointFem3dAttachment (part);
   ax.setFromNodes (nodes, weights);
   mech.addAttachment (ax);

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:

  void setFromNodes (FemNode[] nodes, double[] weights)
  void setFromNodes (Collection<FemNode> nodes, VectorNd weights)
  boolean setFromNodes (Point3d pos, FemNode[] nodes)
  boolean setFromNodes (Point3d pos, Collection<FemNode> nodes)

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

\alpha_{k}=\frac{d_{\text{max}}}{d_{k}+d_{max}} (6.6)

where d_{k} is the distance from node k to pos and d_{\text{max}} 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.

6.4.8 Example: element vs. nodal-based attachments

Figure 6.9: PointFemAttachment loaded into ArtiSynth and run until stable. The top and bottom models are connected to their springs using element and nodal-based attachments, respectively. The nodes associated with each attachment are rendered as white spheres.

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 6.9 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.

1    // Filter to select only elements for which the nodes are entirely on the
2    // positive side of the x-z plane.
3    private class MyFilter extends ElementFilter {
4       public boolean elementIsValid (FemElement e) {
5          for (FemNode n : e.getNodes()) {
6             if (n.getPosition().y < 0) {
7                return false;
8             }
9          }
10          return true;
11       }
12    }
13
14    // Collect and return all the nodes of an FEM model associated with a
15    // set of elements specified by an array of element numbers
16    private HashSet<FemNode3d> collectNodes (FemModel3d fem, int[] elemNums) {
17       HashSet<FemNode3d> nodes = new HashSet<FemNode3d>();
18       for (int i=0; i<elemNums.length; i++) {
19          FemElement3d e = fem.getElements().getByNumber (elemNums[i]);
20          for (FemNode3d n : e.getNodes()) {
21             nodes.add (n);
22          }
23       }
24       return nodes;
25    }
26
27    public void build (String[] args) {
28       MechModel mech = new MechModel ("mech");
29       addModel (mech);
30       mech.setGravity (0, 0, 0); // turn off gravity
31
32       // create and add two FEM beam models centered at the specified locations
33       FemModel3d fem1 = addFem (mech,  0.0, 0.0,  0.25);
34       FemModel3d fem2 = addFem (mech,  0.0, 0.0, -0.25);
35
36       // reconstruct the FEM surface meshes so that they show only elements on
37       // the positive side of the x-y plane. Also, set surface rendering to
38       // show strain values.
39       fem1.createSurfaceMesh (new MyFilter());
40       fem1.setSurfaceRendering (SurfaceRender.Strain);
41       fem2.createSurfaceMesh (new MyFilter());
42       fem2.setSurfaceRendering (SurfaceRender.Strain);
43
44       // create and add the particles for the point-to-point springs
45       // that will apply forces to each FEM.
46       Particle p1 = addParticle (mech,  0.9, 0.0,  0.25);
47       Particle p2 = addParticle (mech,  0.9, 0.0, -0.25);
48       Particle m1 = addParticle (mech,  0.0, 0.0,  0.25);
49       Particle m2 = addParticle (mech,  0.0, 0.0, -0.25);
50
51       // attach spring end-point to fem1 using an element-based marker
52       mech.attachPoint (m1, fem1);
53
54       // attach spring end-point to fem2 using a larger number of nodes, formed
55       // from the node set for elements 22, 31, 40, 49, and 58. This is done by
56       // explicitly creating the attachment and then setting it to use the
57       // specified nodes
58       HashSet<FemNode3d> nodes =
59          collectNodes (fem2, new int[] { 22, 31, 40, 49, 58 });
60
61       PointFem3dAttachment ax = new PointFem3dAttachment (m2);
62       ax.setFromNodes (m2.getPosition(), nodes);
63       mech.addAttachment (ax);
64
65       // finally, create the springs
66       addSpring (mech, /*stiffness=*/10000, p1, m1);
67       addSpring (mech, /*stiffness=*/10000, p2, m2);
68
69       // set the attachments nodes for m1 and m2 to render as white spheres
70       setAttachedNodesWhite (m1);
71       setAttachedNodesWhite (m2);
72       // set render properties for m1
73       RenderProps.setSphericalPoints (m1, 0.015, Color.GREEN);
74    }

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 (6.6).

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 6.9. 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.