6 Finite Element Models

6.10 Fields

In modeling applications, particularly those employing FEM methods, situations often arise where it is necessary to describe numeric quantities that vary over some spatial domain. For example, the stiffness parameters of an FEM material may vary at different points within the volumetric mesh. When modeling muscles, the activation direction vector will also typically vary over the mesh.

ArtiSynth provides field objects which can be used to represent such spatialy varying quantities, together with mechanisms to attach these to properties within FEM materials. Field objects can provide either scalar or vector fields. Scalar fields implement the interface ScalarField and supply the method

   double getValue (Point3d pos)

while vector fields implement VectorField and supply the method

   T getValue (Point3d pos)

where T is any class implementing maspack.matrix.VectorObject. (Most of the fixed-size matrix objects in maspack.matrix, along with MatrixNd and VectorNd, implement VectorObject.) In both cases, the idea is to provide values at arbitrary points over some spatial domain.

Both ScalarField and VectorField are subclassed from Field. The reason for separating scalar and vector fields is simply efficency: having scalar fields work with the primitive type double, instead of the object Double, requires considerably less storage and somewhat less computational effort.

6.10.1 Field types

A field is typically implemented by specifying a finite set of values at discrete locations on an underlying spatial grid and then using different methods to determine values at arbitrary points.

At present, ArtiSynth provides the following different field types:

Grid fields

Implemented by classes that extend either ScalarGridField or VectorGridField, grid fields specify their values at the vertices of a regular 3D grid. Values at an arbitrary position {\bf p} are determined by finding the grid cell containing {\bf p}, and then using trilinear interpolation of the surrounding nodal values.

Nodal fields

Implemented by classes that extend either ScalarNodalField or VectorNodalField, nodal fields specify their values at the nodes of an FEM mesh. Values at an arbitrary position {\bf p} are determined by finding the element containing {\bf p}, and then finding the nodal coordinates (i.e., weights) of {\bf p} within that element and interpolating the value accordingly. Values within a nodal field can also be queried directly either by node, or a set of node numbers and weights:

   T getValue (FemNode3d node)

   T getValue (int[] nodeNums, double[] weights)
Element fields

Implemented by classes that extend either ScalarElementField or VectorElementField, element fields specify their values at the elements of an FEM mesh, with the value assumed to be constant over the element. Values at an arbitrary position {\bf p} are determined by finding the element containing {\bf p} and then returning the value for that element. Values within a nodal field can also be queried directly by element:

   T getValue (FemElement3dBase elem)
Sub-element fields

Implemented by classes that extend either ScalarSubElemField or VectorSubElemField, sub-element fields specify their values at specific points within each element of an FEM mesh, with the points corresponding to the element’s integration points. Values at an arbitrary position {\bf p} are determined by finding the element containing {\bf p}, extrapolating the point values back to the nodes, and then using nodal interpolation. While this is computationally expensive (and it is also possible to create a nodal field that gives the same result), the purpose of sub-element fields is to provide precise values at element integration points, which can be queried directly using the element and point sub-index:

   T getValue (FemElement3dBase elem, int subIdx)

Internally, field features, such as nodes or elements, are generally referred to by their number. This is because node and element numbers are guaranteed to be persistent, in the sense that the number will remain unchanged as long as the node or element is not removed from the FEM model. Element field implementations also need to distingish between volumetric and shell elements, since these are each numbered independently.

Table 6.6 gives a list of all the currently implemented fields. Some of them, such as those prefaced with Vector3d, VectorNd, and MatrixNd, are specialized because they implement special features for their value type, or require specialized construction because the value type does not have a fixed size.

Field component value field type
modelbase.ScalarGridField double scalar grid
modelbase.VectorGridField VectorObject vector grid
femmodels.ScalarNodalField double scalar nodal
femmodels.VectorNodalField VectorObject vector nodal
femmodels.Vector3dNodalField Vector3d vector nodal
femmodels.VectorNdNodalField VectorNd vector nodal
femmodels.MatrixNdNodalField MatrixNd vector nodal
femmodels.ScalarElementField double scalar element
femmodels.VectorElementField VectorObject vector element
femmodels.Vector3dElementField Vector3d vector element
femmodels.VectorNdElementField VectorNd vector element
femmodels.MatrixNdElementField MatrixNd vector element
femmodels.ScalarSubElemField double scalar sub-element
femmodels.VectorSubElemField VectorObject vector sub-element
femmodels.Vector3dSubElemField Vector3d vector sub-element
femmodels.VectorNdSubElemField VectorNd vector sub-element
femmodels.MatrixNdSubElemField MatrixNd vector sub-element
Table 6.6: List of all the field components available in ArtiSynth, with their value and field types.

Each field component also has a default value, which is returned for any query for which the specified position lies outside the field’s domain. For a grid field, this means lying outside the bounds of the grid. For an FEM-based field, this means lying outside the FEM mesh. The default value can be specified in the field’s constructor. If not specified, it is usually set to 0.

Having a default value also means that explicit values do not have to specified at all of the features (e.g., nodes or elements) used to define the field. If an explcit value is missing, the default value is used instead.

6.10.2 Adding fields to a model

Fields are added to an ArtiSynth model by first creating an instance of the desired field class, populating it with values (at vertices, nodes, or elements, as required by the field), and then adding it to the model. Grid fields can generally be added anywhere, whereas FEM-based fields should be added to the fields sublist of their associated FEM model (grid fields can be added there as well).

As a simple example, one can construct a scalar nodal field, associated with an FEM model named fem, like this:

   FemModel3d fem;
   // ... build the fem ...
   ScalarNodalField field = new ScalarNodalField ("stiffness", fem, 0);
   for (FemNode3d n : fem.getNodes()) {
      double value = ... // compute value for the node
      field.setValue (n, value);
   }
   fem.addField (field); // add the field to the fem

Each field type has several constructors available. The one used above takes a component name ("stiffness"), the FEM modal associated with the field, and a default value (0 in this case). When setting field values, it is often easiest to simply specify the node (although the node’s number (i.e., node.getNumber()) can be specified as well. After being created and initialized, the field is added to the FEM model using its addField() method.

As noted in the previous section, it is not necessary to specify values for all the features of a field (e.g., nodes for a nodal field, or elements for an element field). For features with undefined values, the default value will be used instead.

As another example, consider constructing a vector element field, with two-dimensional values characterized by Vector2d. An example of this is as follows:

   VectorElementField<Vector2d> field =
      new VectorElementField<> ("uv", Vector2d.class, fem);
   for (FemElement3d e : fem.getElements()) {
      Vector2d value = ... // compute value for the element
      field.setValue (e, value);
   }
   fem.addField (field); // add the field to the fem

General vector fields are declared using type parameterization, in which the class type of the vector must be appended to the declaration between angle brackets (< >). Also, the constructor also requires the class type to be specified as an argument. In the example above, the field is constructed with a name ("uv"), value class type (Vector2d.class), and the FEM model. Since no default value is specified, a value of 0 will be assumed.

Since it is widely used, the vector type Vector3d has its own predefined field classes, which don’t require parameterization and may also contain special methods and features relavant for 3D vectors. A Vector3d nodal field can be declared and created as follows:

    Vector3dNodalField field = new Vector3dNormalField ("normals", fem);

Since the vector type is explicitly “wired in”, no type parameterization is required, and Vector3d.class does not need to be specified to the constructor.

The vector types VectorNd and MatrixNd also have special predefined classes, because these components do not have a predetermined size and so it is necessary to specify this size information when constructing the field. Some constructors for fields with these types look like:

   MatrixNdElementField field =
      new MatrixNdElementField ("matrixField", 3, 6, fem);
   VectorNdNodalField field =
      new VectorNdNodalField ("params", 7, fem);

Both of these contain size information after the name argument: 3\times 6 for the matrix field, and 7 for the vector field. Again, no default value is specified, so 0 will be assumed.

Lastly, we consider constucting a sub-element field. This is normally done when we need precisely computed information at each of the element’s integration points, and we can’t assume that nodal interpolation will give an accurate enough approximation. A general template for constructing a sub-element field might look like this:

   ScalarSubElemField field = new ScalarSubElemField ("precise", fem, -1);
   for (FemElement3d e : fem.getElements()) {
      IntegrationPoint3d[] ipnts = e.getAllIntegrationPoints();
      for (int k=0; k<ipnts.length; k++) {
         double value = ... // compute value at integration point k
         field.setValue (e, k, value);
      }
   }
   fem.addField (field); // add the field to the fem

Here, we create a scalar sub-element field named "precise" with a default value of -1. Then we iterate first through the FEM elements, and then through each element’s integration points, computing values appropriate to each one. Since we are likely to need information from the integration points themselves to compute the value, we use e.getAllIntegrationPoints() to obtain a list of them.

Assigning values to a sub-element field, it is important to use the element’s getAllIntegrationPoints() method instead of getIntegrationPoints(), since the former returns a list of all integration points, including the special warping point which is contained in the element’s center and is used by the solver under some circumstances. Likewise, the number of sub-element values required is given by numAllIntegrationPoints().

When computing sub-element field values, if it is necessary to obtain either the rest or current (spatial) position for an integration point, these can be obtained as follows:

      IntegrationPoint3d ipnt;
      Point3d pos = new Point3d();
      FemElement3d elem;
      ...
      // compute spatial position:
      ipnt.computePosition (pos, elem.getNodes());
      // compute rest position:
      ipnt.computeRestPosition (pos, elem.getNodes());

6.10.3 Binding to material propertues

Once a field has been created and added to the model, one may bind it to certain properties of a FemMaterial. When this is done, then whenever those properties are queried internally by the solver at the integration points for each FEM element, the field value at that integration point is used instead of the regular property value.

As long as a property XXX is bound to a field, its regular value will still appear (and can be set) through widgets in control or property panels, or via the getXXX() and setXXX() accessors. However, the regular value won’t be used internally by the FEM simulation.

To bind a property to a field, it is necessary that

  1. 1.

    The type of the field matches the value of the property;

  2. 2.

    The property is itself bindable.

If the property has a double value, then it can be bound to any ScalarField. Otherwise, if the property has a value T, where T is an instance of maspack.matrix.VectorObject, then it can be bound to any VectorField<T>.

Bindable properties export two methods with the signature

   setXXXField (T field, boolean useRestPos)
   T getXXXField ()

where XXX is the property name and T is an appropriate field type. The second arugment, useRestPos, is relevant only for grid fields, and indicates whether the grid should be queried using the integration point’s rest position, or current spatial position. For FEM-based fields, this distinction is irrelevant, since the field is assumed to always move with the FEM.

For example, consider the YoungsModulus property of a NeoHookeanMaterial material. This has a double value, and is bindable, and so can be bound to a ScalarField as follows:

   FemModel3d fem;
   ScalarNodalField field;
   NeoHookeanMaterial mat;
   ... other initialization ...
   mat.setYoungsModulusField (field, false);  // bind to the field
   fem.setMaterial (mat);   // set the material in the FEM model

It is important to perform field bindings on materials before they are set in an FEM model (or one of its subcomponents, such as MuscleBundles). That’s because the setMaterial() method copies the input material, and so any settings made on it afterwards won’t be seen by the FEM:

   // works: field will be seen by the copied version of ’mat’
   mat.setYoungsModulusField (field, false);
   fem.setMaterial (mat);

   // does NOT work: field not seen by the copied version of ’mat’
   fem.setMaterial (mat);
   mat.setYoungsModulusField (field, false);

To unbind YoungsModulus, one can do

   mat.setYoungsModulusField (null, false);

where the value of the useRestPos argument is ignored.

6.10.4 Example: FEM with variable stiffness

Figure 6.15: VariableStiffness model after being run in ArtiSynth.

A simple model demonstrating a stiffness that varies over an FEM mesh is defined in

  artisynth.demos.tutorial.VariableStiffness

It consists of a simple thin hexahedral beam with a linear material for which the Youngs modulus E is made to vary nonlinearly along the x axis of the rest position according to the formula

E=\frac{10^{8}}{1+1000\,x^{3}} (6.10)

The model’s build method is given below:

1    public void build (String[] args) {
2       MechModel mech = new MechModel ("mech");
3       addModel (mech);
4
5       // create regular hex grid FEM model
6       FemModel3d fem = FemFactory.createHexGrid (
7          null, 1.0, 0.25, 0.05, 20, 5, 1);
8       fem.transformGeometry (new RigidTransform3d (0.5, 0, 0)); // shift right
9       fem.setDensity (1000.0);
10       mech.addModel (fem);
11
12       // fix the left-most nodes
13       double EPS = 1e-8;
14       for (FemNode3d n : fem.getNodes()) {
15          if (n.getPosition().x < EPS) {
16             n.setDynamic (false);
17          }
18       }
19       // create a scalar nodel field to make the stiffness vary
20       // nonlinearly along the rest position x axis
21       ScalarNodalField stiffnessField = new ScalarNodalField(fem, 0);
22       for (FemNode3d n : fem.getNodes()) {
23          double s = 10*(n.getRestPosition().x);
24          double E = 100000000*(1/(1+s*s*s));
25          stiffnessField.setValue (n, E);
26       }
27       fem.addField (stiffnessField);
28       // create a linear material, bind its Youngs modulus property to the
29       // field, and set the material in the FEM model
30       LinearMaterial linearMat = new LinearMaterial (100000, 0.49);
31       linearMat.setYoungsModulusField (stiffnessField, /*useRestPos=*/true);
32       fem.setMaterial (linearMat);
33
34       // set some render properties for the FEM model
35       fem.setSurfaceRendering (SurfaceRender.Shaded);
36       RenderProps.setFaceColor (fem, Color.CYAN);
37    }

Lines 6-10 create the hex FEM model and shift it so that the left side is aligned with the orihin, while lines 12-17 fix the leftmost nodes. Lines 21-27 create a scalar nodal field for the Youngs modulus, with lines 23-24 computing E according to (6.10). The field is then bound to a linear material which is then set in the model (lined 30-32).

The example can be run in ArtiSynth by selecting All demos > tutorial > VariableStiffness from the Models menu. When it is run, the beam will bend under gravity, but mostly on the right side, due to the much higher stiffness on the left (Figure 6.15).

6.10.5 Example: specifying FEM muscle directions

Figure 6.16: (Top) RadialMuscle model after being loaded into ArtiSynth and run with its excitation set to 0. (Bottom) RadialMuscle with its excitation set to around .375.

Another example involves using a Vector3d field to specify the muscle activation directions over an FEM model and is defined in

  artisynth.demos.tutorial.RadialMuscle

When muscles are added using MuscleBundles ( Section 6.9), the muscle directions are stored and handled internally by the muscle bundle itself. However, it is possible to add a MuscleMaterial directly to the elements of an FEM model, using a MaterialBundle (Section 6.8), in which case the directions need to be set explicitly using a field.

The model’s build method is given below:

1    public void build (String[] args) {
2       MechModel mech = new MechModel ("mech");
3       addModel (mech);
4
5       // create a thin cylindrical FEM model with two layers along z
6       double radius = 0.8;
7       FemMuscleModel fem = new FemMuscleModel ("radialMuscle");
8       mech.addModel (fem);
9       fem.setDensity (1000);
10       FemFactory.createCylinder (fem, radius/8, radius, 20, 2, 8);
11       fem.setMaterial (new NeoHookeanMaterial (200000.0, 0.33));
12       // fix the nodes close to the center
13       for (FemNode3d node : fem.getNodes()) {
14          Point3d pos = node.getPosition();
15          double radialDist = Math.sqrt (pos.x*pos.x + pos.y*pos.y);
16          if (radialDist < radius/2) {
17             node.setDynamic (false);
18          }
19       }
20       // compute a direction field, with the directions arranged radially
21       Vector3d dir = new Vector3d();
22       Vector3dElementField dirField = new Vector3dElementField (fem);
23       for (FemElement3d elem : fem.getElements()) {
24          elem.computeCentroid (dir);
25          // set directions only for the upper layer elements
26          if (dir.z > 0) {
27             dir.z = 0; // remove z component from direction
28             dir.normalize();
29             dirField.setValue (elem, dir);
30          }
31       }
32       fem.addField (dirField);
33       // add a muscle material, and use it to hold a simple force
34       // muscle whose ’restDir’ property is attached to the field
35       MaterialBundle bun = new MaterialBundle ("bundle",/*all elements=*/true);
36       fem.addMaterialBundle (bun);
37       SimpleForceMuscle muscleMat = new SimpleForceMuscle (500000);
38       muscleMat.setRestDirField (dirField, true);
39       bun.setMaterial (muscleMat);
40
41       // add a control panel to control the excitation
42       ControlPanel panel = new ControlPanel();
43       panel.addWidget (bun, "material.excitation", 0, 1);
44       addControlPanel (panel);
45
46       // set some rendering properties
47       fem.setSurfaceRendering (SurfaceRender.Shaded);
48       RenderProps.setFaceColor (fem, new Color (0.6f, 0.6f, 1f));
49    }

Lines 5-19 create a thin cylindrical FEM model, centered on the origin, with radius r and height r/8, consisting of hexes with wedges at the center, with two layers of elements along the z axis (which is parallel to the cylinder axis). Its base material is set to a neo-hookean material. To keep the model from falling under gravity, all nodes whose distance to the z axis is less than r/2 are fixed.

Next, a Vector3d field is created to specify the directions, on a per-element basis, for the muscle material which will be added subsequently (lines 21-32). While we could create an instance of VectorElementField<Vector3d>, we use Vector3dElementField, since this is available and may provide additional functionality (such as the ability to render the directions). Directions are set to lie outward in a radial direction perpendicular to the z axis, and since the model is centered on the origin, they can be computed easily by first computing the element centroids, removing the z component, and then normalizing. In order to give the muscle action an upward bias, we only set directions for elements in the upper layer. Direction values for elements in the lower layer will then automatically have a default value of 0, which will cause the muscle material to not apply any stress.

We next add to the model a muscle material whose directions will be determined by the field. To hold the material, we first create and add a MaterialBundle which is set to act on all elements (line 35-36). Then we set this bundle’s material to SimpleForceMuscle, which simply adds a stress along the muscle direction that equals the excitation value times the value of its maxStress property, and bind the material’s restDir property to the direction field (lines 37-39).

Finally, we create and add a control panel to allow interactive control over the muscle material’s excitation property (lines 42-44), and set some rendering properties for the FEM model.

The example can be run in ArtiSynth by selecting All demos > tutorial > RadialMuscle from the Models menu. When it is first runs, it falls around the edges under gravity (Figure 6.16, top). Applying an excitation causes a radial contraction which pulls the edges upward and, if high enough, causes then to buckle (Figure 6.16, bottom).