Controlling Anisotropy
in Mass-Spring Systems
David Bourguignon and Marie-Paule Cani
iMAGIS-GRAVIR/IMAG-INRIA
iMAGIS is a joint research project of CNRS/INRIA/UJF/INPG
E-mail:
f
David.Bourguignon
j
Marie-Paule.Cani
g
@imag.fr
http://www-imagis.imag.fr/
Abstract. This paper presents a deformable model that offers control of the
isotropy or anisotropy of elastic material, independently of the way the object is
tiled into volume elements. The new model is as easy to implement and almost
as efficient as mass-spring systems, from which it is derived. In addition to con-
trolled anisotropy, it contrasts with those systems in its ability to model constant
volume deformations. We illustrate the new model by animating objects tiled
with tetrahedral and hexahedral meshes.
1 Introduction
Mass-spring systems have been extensively used in Computer Graphics over the last fif-
teen years, and are still very popular. Easier to implement and faster than finite element
methods, these systems allow animation of dynamic behaviors. They have been ap-
plied to the animation of inanimate bodies such as cloth or soft material [17, 10, 15, 5]
and to the animation of organic active bodies such as muscles in character anima-
tion [11, 1, 12].
Both isotropic and anisotropic elastic materials can be found among the objects
to animate. For instance, a rubber toy is isotropic, while most natural objects (animal
organs, plants) are strongly anisotropic, due to their fiber structure and/or the composite
materials they are made of. One of the main drawbacks of mass-spring systems is
that neither isotropic nor anisotropic materials can be generated and controlled easily.
Another problem is that most of the materials found in nature maintain a constant or
quasi-constant volume during deformations (this is well known for muscles, but also
holds for inanimate materials). Mass-spring models do not have this property.
1.1 Background
Animating an elastic object using a mass-spring system usually consists of discretizing
the object with a given 3D mesh, setting point masses on the mesh nodes and damped
springs on the mesh edges. Then, most implementations simply integrate point dynam-
ics equations for each mass from the set of applied forces due to the mesh deformation
at the previous time step [10].
Well-known advantages of mass-spring systems include their ability to generate
dynamic behaviors, while finite elements methods are generally used in the static case,
and their ability to handle both large displacements and large deformations.
Among the intrinsic limitations of mass-spring systems, one of the main problems is
parameter setting. Computing the masses in order to set up a homogeneousmaterial can
be done by computing each mass according to the volume of the Voronoi region around
it [4]. However, there is no easy solution for spring parameters. Since damped springs
are positioned along the edges of a given volume mesh, the geometrical and topological
structure of this mesh strongly influences the material behavior. A consequence of this
problem is that changing the mesh density during the simulation while maintaining the
same global mechanical properties is very difficult [7].
If all springs are set to the same stiffness, the mesh geometry may generate unde-
sired anisotropy, as shown in Fig. 1.a. The undesired behavior disappears when hexa-
hedral elements aligned with the forces directions are used (Fig. 1.b). Of course, if the
tiling of the object volume was computed from the triangulation of random uniformly-
distributed sample points, the unwanted anisotropy problem would tend to disappear
when the density of the mesh increases. However, using an extremely dense mesh
would reduce efficiency.
Approximating a desired behavior using a given mesh can be achieved, as in [4, 9],
by using optimization to tune individual spring stiffnesses. This technique could be
used, in theory, for generating both isotropic and anisotropic behaviors. However, due
to the large computational cost, this method has only been tested in the 2D case [4].
The most common approach to control the behavior of a mass-spring system, at
least along a few “directions of interest”, is to specifically design the mesh in order
to align springs on these specific directions, such as in Fig. 1.b. This was done for
instance in Miller’s “snakes and worms” models [11] and in the muscle model of Ng and
Fiume [12], where some of the springs were aligned with the muscle fibers and the rest
were set perpendicular to them. Unfortunately, manually creating such meshes would
be time consuming in the general case, where fiber directions generating anisotropy
vary in an arbitrary way inside the object. We are rather looking for an approach that
uses a 3D mesh obtained, for example, with a commercial meshing package (such as
GHS3D [16]) fed with a 3D surface mesh, and still displays the deformable model
behavior, with specified properties in specific directions.
00
00
00
11
11
11
shear spring
point mass
structural spring
000
000
000
000
111
111
111
111
abc d
Fig. 1. Mass-spring systems drawbacks. At left, comparison between two meshes undergoing a
downward pull at their bottom end while their top end is fixed. We observe undesired anisotropy
in the tetrahedral mass-spring system (a), but not in the hexahedral mesh with springs aligned
in the gravity and pull force directions (b). At right, equilibrium state of a cantilever beam,
which left end is fixed, under force of gravity (c). All things being equal, the mass-spring system
considered (tetrahedral mesh) is unable to sustain flexion, as opposed to our model (Fig. 5.c).
The spring configurations used for tetrahedral and hexahedral meshes are given in (d).
1.2 Overview
This paper presents an alternative model to classical mass-spring systems that enables
one to specify isotropic or anisotropic properties of an elastic material, independently
from the 3D mesh used for sampling the object. The approach we use is still related
to mass-spring systems, in the sense that we animate point masses subject to applied
forces. However, the forces acting on each mass are derived from the anisotropic be-
havior specified for each of the volume elements that are adjacent to it.
Since there are no springs along the mesh edges, the geometry and topology of the
mesh do not restrict the simulated behavior. Moreover, constant volume deformations
can be obtained easily, by adding extra forces. We illustrate this on both tetrahedral and
hexahedral meshes. Our results show that computation time remains low, while more
controllable behaviors are achieved.
2 Modeling Anisotropy
Our aim is to specify the mechanical properties of the material independently from
the mesh geometry and topology. In usual mass-spring systems, internal forces acting
inside the material are approximated exclusively by forces acting along the edges of the
mesh (i.e. along the springs). This is the reason for the undesired anisotropy problem
described earlier, and for the difficulty in specifying desired anisotropic properties.
The basic idea of our method is to let the user define, everywhere in the object,
mechanical characteristics of the material along a given number of axes corresponding
to orientations of interest at each current location. All internal forces will be acting
along these axes instead of acting along the mesh edges. For instance, in the case of
organic materials such as muscles, one of the axes of interest should always correspond
to the local fiber orientation.
Since the object is tiled using a mesh, axes of interest and the associated mechanical
properties are specified at the barycenter of each volume element inside the mesh. We
currently use three orthogonal axes of interest. The possible use of a larger number of
axes will be discussed in Section 6.
2.1 General Scheme
During deformations of the material, the three axes of interest, of given initial orien-
tation, evolve with the volume element to which they belong. In order to be able to
know their position at each instant, we express the position of the intersection point of
one axis with one of the element faces as a linear combination of the positions of the
vertices defining the face. The corresponding interpolation coefficients are computed
for each face in the rest position (see Figures 2 and 4).
Given the position of the point masses of a volume element, we are thus able to
determine the coordinates of the six intersection points and consequently the three axes
that constitutes the local frame, up to the precision of our linear interpolation.
From the deformation of the local frame, we can deduce the resulting forces on
each intersection point. Then, for a given face, we can compute the force value on each
point mass belonging to this face by “inverse” interpolation of the force value at the
intersection point. The interpolation coefficients previously defined are therefore also
considered as weighting coefficients of the force on each point mass.
2.2 Forces Calculations
Damped springs with associated stiffness and damping coefficients are used to model
stretching characteristics along each axis of interest. In order to specify shearing prop-
erties, angular springs are added between each pair of axes. Rest lengths and rest angles
are pre-computed from the initial position of the object that defines its rest shape. The
equations we use for these springs are detailed below.
Axial damped spring. The spring forces f
1
and f
2
between a pair of intersection points
1 and 2 at positions x
1
and x
2
with velocities v
1
and v
2
are
f
1
=
;
k
s
(
k
l
21
k;
r
)+
k
d
˙
l
21
l
21
k
l
21
k
l
21
k
l
21
k
;
f
2
=
;
f
1
;
where l
21
=
x
1
;
x
2
, r is the rest length,
˙
l
21
=
v
1
;
v
2
is the time derivative of l
21
, ks
and kd are respectively the stiffness and damping constants.
Angular spring. The spring forces
(
f
1
;
f
2
)
and
(
f
3
;
f
4
)
between two pairs of intersec-
tion points
(
1
;
2
)
and
(
3
;
4
)
are
f
1
=
;
k
s
l
21
l
43
k
l
21
kk
l
43
k
;
c
l
43
k
l
43
k
;
f
2
=
;
f
1
;
f
3
=
;
k
s
l
21
l
43
k
l
21
kk
l
43
k
;
c
l
21
k
l
21
k
;
f
4
=
;
f
3
;
where l
21
=
x
1
;
x
2
and l
43
=
x
3
;
x
4
, c is the cosine of the rest angle between l
21
and
l
43
, ks is the stiffness constant.
Here, two approximations are made: first, we assume a small variation of the angle
and take the variation of the angle’s cosine instead; second, we consider it sufficient
to use as unit vector the other vector of the pair, instead of a vector normal to the one
considered, in the plane where the angle is measured. These two approximations gave
good results in practice. Furthermore, we found no necessity to use damped angular
springs.
3 Application to Tetrahedral Meshes
Many objects in Computer Graphics are modeled using triangular surface meshes. Gen-
erating a 3D mesh from such a description, using tools like GHS3D [16] yields to tetra-
hedral volume meshes. This section details our method in this case.
Fig. 2 depicts a tetrahedral element, with the associated frame defining the three
axes of interest. We express the position x
P
of point P as a function of the positions of
vertices A, B and C of the given face, using barycentric coordinates:
x
P
=
α x
A
+
β x
B
+
γ x
C
(e.g. if α
=
1 and β
=
γ
=
0, we get x
P
=
x
A
). Therefore, a force f
P
applied to point P
is split into forces αf
P
, βf
P
and γf
P
, respectively applied on points A, B and C.
We can note that since the elementary volume has four faces, and since there are
three axes of interest defining six intersection points, two such points may lie on the
same face of the volume. This has not been problematic in practice, since forces applied
on mesh nodes are correctly weighted.
P
α
BA
C
γ
BA
C
β
Fig. 2. Tetrahedral element. A point mass is located at each vertex. A local frame is defined
at the barycenter of the element (left). Each axis is characterized by the barycentric coordinates
α, β and γ (with α
+
β
+
γ
=
1) of its two intersection points (right, for a given face). These
coordinates are easily obtained using an area ratio.
3.1 Volume Preservation
Animating constant volume deformations with a classical mass-spring system is not
straightforward. For these systems, forces are only applied along the edges of each
volume element, while maintaining a constant volume basically requires adding radial
forces or displacements, as shown by Promayon et al. [14].
To simply ensure volume preservation, we propose a volume force formulation
adapted to tetrahedral volume element. It is loosely related to soft volume-preservation
constraint of Lee et al. [8].
Let us define x
B
the position of the barycenter of the tetrahedral element, with
x
B
=
1
4
3
∑
i
=
0
x
i
where x
i
is the position of the ith vertex. Then, we define the force applied on the jth
vertex as
f
j
=
;
k
s
"
3
∑
i
=
0
k
x
i
;
x
B
k;
3
∑
i
=
0
k
x
i
;
x
B
k
t
=
0
#
x
j
;
x
B
k
x
j
;
x
B
k
where k
s
is the constraint stiffness and
∑
3
i
=
0
k
x
i
;
x
B
k
t
=
0
is the rest length of this “vol-
ume spring”. It was not necessary to add damping forces with this constraint.
This method gave satisfactory results in pratice, since we get less than 1
:
5% volume
variation in our experiment (see Fig. 3), but results depend on the material parameters
chosen and the type of experiment conducted. In applications where these volume
variations are considered too high, volume preservation could be enforced directly as a
hard constraint like in Witkin’s work [20, 19].
4 Application to Hexahedral Meshes
The use of hexahedral meshes is not as common as tetrahedral ones, since the geometry
they can define is more limited. However, these meshes may be useful for animating
objects modeled using voxels [2]. This kind of data, with information about material
characteristics specified in each voxel (possibly including anisotropy), may be provided
by medical imaging applications.
Applying the general method presented in Section 2.1 to hexahedral meshes is
straightforward. Fig. 4 depicts an hexahedral element, with the associated frame defin-
ing the three axes of interest. We express the position x
P
of point P as a function of
