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 spatially 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
while vector fields implement VectorField and supply the method
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 efficiency: having scalar fields work with the primitive type double, instead of the object Double, requires considerably less storage and somewhat less computational effort.
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:
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 are determined by finding the grid cell containing , and then using trilinear interpolation of the surrounding nodal values.
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 are determined by finding the element containing , and then finding the nodal coordinates (i.e., weights) of 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)
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 are determined by finding the element containing and then returning the value for that element. Values within a nodal field can also be queried directly by element:
T getValue (FemElement3dBase elem)
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 are determined by finding the element containing , 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 distinguish 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|
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 explicit value is missing, the default value is used instead.
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:
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:
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 relevant for 3D vectors. A Vector3d nodal field can be declared and created as follows:
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:
Both of these contain size information after the name argument: for the matrix field, and 7 for the vector field. Again, no default value is specified, so 0 will be assumed.
Lastly, we consider constructing 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:
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:
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
The type of the field matches the value of the property;
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
where XXX is the property name and T is an appropriate field type. The second argument, 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:
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
where the value of the useRestPos argument is ignored.
A simple model demonstrating a stiffness that varies over an FEM mesh is defined in
It consists of a simple thin hexahedral beam with a linear material for which the Young’s modulus is made to vary nonlinearly along the x axis of the rest position according to the formula
The model’s build method is given below:
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 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).
Another example involves using a Vector3d field to specify the muscle activation directions over an FEM model and is defined in
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:
Lines 5-19 create a thin cylindrical FEM model, centered on the origin, with radius and height , consisting of hexes with wedges at the center, with two layers of elements along the 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 axis is less than 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 axis, and since the model is centered on the origin, they can be computed easily by first computing the element centroids, removing the 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).