7 Fields

7.5 Binding material properties

A principal purpose of field components is to enable certain FEM material properties to vary spatially over the FEM geometry. Many material properties can be bound to a field, so that when they 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. Likewise, some properties of contact force behaviors, such as the YoungsModulus and thickness properties of LinearElasticContact (Section 8.7.3) can be bound to a field that varies over the contact mesh geometry.

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 ScalarFieldComponent. Otherwise, if the property has a value T, where T is an instance of maspack.matrix.VectorObject, then it can be bound to a vector field.

Bindable properties export two methods with the signature

⬇
   setXXXField (F field)
   F getXXXField ()

where XXX is the property name and F is an appropriate field type.

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 its getXXX() and setXXX() accessors. However, the regular value won’t be used internally by the FEM simulation.

For example, consider the YoungsModulus property, which is present in several FEM materials, including LinearMaterial and NeoHookeanMaterial. This has a double value, and is bindable, and so can be bound to a ScalarFieldComponent as follows:

⬇
   FemModel3d fem;
   ScalarNodalField field;
   NeoHookeanMaterial mat;
   ... other initialization ...
   mat.setYoungsModulusField (field);  // 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);
   fem.setMaterial (mat);

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

To unbind YoungsModulus, one can call

⬇
   mat.setYoungsModulusField (null);

Usually, FEM fields are used to bind properties in FEM materials while mesh fields are used for contact properties, but other fields can be used, particularly grid fields. To understand this a bit better, we discuss briefly some of the internals of how binding works. When evaluating stresses, FEM materials have access to a FemFieldPoint object that provides information about the integration point, including its element, nodal weights, and integration point index. This can be used to rapidly evaluate values within nodal, element and sub-element FEM fields. Likewise, when evaluating contact forces, contact materials like LinearElasticContact have access to a MeshFieldPoint object that describes the contact point position as a weighted sum of mesh vertex positions, which can then be used to rapidly evaluate evaluate values within vertex or face mesh fields. However, all fields, regardless of type, implement methods for determining values at both FemFieldPoints and MeshFieldPoints:

⬇
  T getValue (FemFieldPoint fp)
  T getValue (MeshFieldPoint mp)

If necessary, these methods simply fall back on the getValue(Point3d) method, which is defined for all fields and works relatively fast for grid fields.

There are some additional details to consider when binding FEM material properties:

  1. 1.

    When binding to a grid field, one has the choice of whether to use the grid to represent values with respect to the FEM’s rest position or spatial position. This can be controlled by setting the useFemRestPositions property of the grid field, the default value for which is true.

  2. 2.

    When binding the bulkModulus property of an incompressible material, it is best to use a nodal field if the FEM model is using nodal-based soft incompressibility (i.e., the model’s softIncompMethod property is set to either NODAL or AUTO; Section 6.7.3). That’s because soft nodal incompressibility require evaluating the bulk modulus property at nodes instead of integration points, which is much easier to do with a nodal-based field.

7.5.1 Example: FEM with variable stiffness

Figure 7.1: 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 Young’s 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}} (7.0)

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 origin, while lines 12-17 fix the leftmost nodes. Lines 21-27 create a scalar nodal field for the Young’s modulus, with lines 23-24 computing E according to (7.0). The field is then bound to a linear material which is then set in the model (lines 30-32).

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

7.5.2 Example: specifying FEM muscle directions
Figure 7.2: (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);
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 provides 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 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 run, it falls around the edges under gravity (Figure 7.2, top). Applying an excitation causes a radial contraction which pulls the edges upward and, if high enough, causes then to buckle (Figure 7.2, bottom).

Bibliography