Sofa Mechanics
Sofa Mechanics
Sofa Mechanics
Abstract Mesh-based techniques are well studied and established methods for
solving continuum biomechanics problems. When the problem at hand involves
extreme deformations or artificial discontinuities, meshless methods provide several advantages over the mesh-based methods. This work discusses the Moving
Least Square approximation-based meshless collocation method for simulating deformable objects and presents a verification technique that is based on the
Hertzian theory of non-adhesive elastic contact. The effectiveness of the Hertzian
contact theory as a means for verification was first tested and proven through a
well-established FEM code, FEBio. The meshless method was implemented as a
reusable component for the SOFA framework, an open source software library for
real-time simulations. Through experimentation, the Hertzian theory has been
tested against SOFA hexahedral FEM and the meshless models within the SOFA
framework. Convergence studies and L2 error curves are provided for both models. Experimental results demonstrated the effectiveness of the implementation of
the meshless method.
1 Introduction
Soft tissue models have a wide range of application areas with particular focus
on real-time medical simulations. As a result, accurate interactive modeling of soft
tissue is an important and well-established research field in continuum biomechanics studies.
A continuum model typically relies on an underlying mesh structure either in
2D or 3D depending on the nature and the requirements of the problem. A
breadth-first classification of mesh-based continuum models includes mass-spring
networks [1], Finite Element Methods [2], Finite Volume Methods [3], and Finite
Difference Methods [4].
Originally developed for molecular dynamics problems and astrophysical simulations, meshless (mesh-free) methods offer an appealing alternative to meshbased methods when the problem involves large deformations and imposed dis-
continuities such as cracks and cuts. In this work, the meshless method is implemented as a component under the open-source SOFA library, which focuses on
real-time interactive simulations with an emphasis on medical simulations. In conjunction with the interactive simulation requirement, we assume the linear property of the soft tissue model within a specific range of the strain-stress curve. This
assumption is parallel with the small strain and linear material assumptions of the
Hertzian contact theory.
There have been numerous approaches to modeling mechanics of soft tissues.
Among those is the early work of Sederberg and Parry [5] that presented a technique for deforming solid geometric models in an intuitive free-form manner. The
deformations were based on interpolating trivariate Bernstein polynomials, and
could be applied either globally or locally with volume preservation. Free-form
deformation is an approximate and simple method for deforming solid objects;
however the lack of physical basis is grounds for excluding it as an option for realistic medical simulations. The work of Frisken-Gibson [6] modified the traditional
voxel based representation of volumetric objects and presented a linked volume
representation that was capable of handling interactive object manipulations such
as carving, cutting, tearing, and joining, but still unable to produce physically realistic results. An alternative to volumetric methods is to use mass-spring models [7]
and membrane based approximations that utilize spring meshes [8]. A spring mesh
is composed of vertices and edges, in which each edge is realized as a spring that
connects vertices pair-wise, and each vertex is idealized as a point mass. Although
spring meshes employ physical equations like Hookes law, it is difficult to reproduce specific elastic material properties even with very careful distribution of
spring stiffness through the mesh and suffer from poor numerical stability.
The early work of Bro-Nielsen discussed a fast adaptation of finite element
modeling to satisfy speed and robustness requirements in a surgical simulation setting [9]. In this work, the body part was modeled as a 3D linear elastic solid that
consisted of particles, which were deformed into a new shape when forces were
applied to the elastic solid. The author incorporated a technique called condensation in order to achieve interactive simulation. In the finite element modeling context, condensation translates into obtaining a more compact version of the system
model by rearranging the terms of the matrix equations. In this work, the author
condensed the equations by only considering the finite element nodes on the surface of the model.
A number of recent techniques have addressed the fidelity versus efficiency
trade-off. Another important work in the area is the finite element model based on
Total Lagrangian Explicit Dynamics (TLED) by Miller et al. [10]. The difference
between the TLED based finite element model and other approaches is using the
original reference configuration of the object to calculate the stress and strain tensors during a simulation step. As a result of expressing computations in the reference coordinates, the authors were able to pre-compute spatial derivatives. The
pre-computation of the spatial derivatives leads to efficiency in terms of computational resources, while being capable of handling geometric and material non-
linearities. In their work, the authors employed the central difference based explicit integration rather than the implicit integration scheme. With this choice, they
were able to avoid solving the set of non-linear algebraic equations that are required by the implicit integration at each time step. However, the use of explicit
integration brings limitation on the time step size. The authors justified their implementation choice by stating that the relatively lower stiffness (Youngs modulus) value of the soft tissue relaxes the time step limitation considerably compared
to the typical simulations involving more stiff material like steel or concrete.
Another attempt to increase the computational efficiency of the elastic model in
the context of interactive simulation was discussed in the work of Marchesseau et
al. [11]. The authors presented a new discretization method called Multiplicative
Jacobian Energy Decomposition (MJED), which allows the simulation to assemble the stiffness matrix of the system faster than the traditional Galerkin FEM
formulation. The authors reported up to five times faster computations for the St.
Venant Kirchoff materials. For validation purposes, the authors compared the
MJED approach to the traditional FEM implementation in the SOFA framework
[12], referred to as the standard FEM implementation.
FEM techniques have dominated the field of computational mechanics in the
past several decades. They have been widely used for modeling physical phenomena such as elasticity, heat transfer, and electromagnetism and they heavily rely on
the assumption of a continuous domain. When the problem domain no longer
complies with this continuum assumption, the rationale behind using an FEM
based solution disappears. FEM is also not well suited to problems involving extreme mesh distortions that result in degenerate element shapes, moving discontinuities that do not align with the element edges such as propagating cracks, and
advanced material transformations such as melting of a solid or freezing. To address these issues, significant interest has been developed towards a different class
of methods for solving PDEs, namely meshless or mesh-free methods [13, 14].
The very first meshless method dated back to 1977 [15] and proposed a smoothed
particle hydrodynamics (SPH) method that was used to model theoretical astrophysical phenomena such as galaxy formation, star formation, stellar collisions,
and dust clouds. Its meshless Lagrangian nature allowed diverse usage areas besides astrophysics such as fluid flow, ballistics, volcanology, and oceanography
[16].
Although the SPH method eliminates the necessity of a mesh structure and allows the solution of unbounded problems, it also has its limitations. Because of its
simplistic kernel based approximation scheme, it fails to reproduce even first order
polynomials, resulting in severe consistency problems [13]. To alleviate this problem, methods that utilize moving least squares (MLS) approximations have been
developed. The first work that used MLS approximations in a Galerkin method is
the work of Nayroles et al. [17], which was refined by Belytschko et al. [18] and
named Element-Free Galerkin (EFG) method. This class of methods, different
from the SPH method, use shape functions in approximations that are essentially
corrected versions of compact supported weight functions. The shape functions
2 Algorithm
with rij
x j xi
hi
(1)
ing and central particles, respectively and hi is the support radius of the central
particle i.
The mass and density of a meshless node are assigned at the beginning and
kept fixed throughout the simulation. The mass values are initialized with
(2)
mi sri 3 ,
th
where is the material density value, ri is the average distance of the i node to
its k-nearest neighbors, and s is a scaling factor that is chosen so that the average
of the assigned densities is close to the actual material density. The assigned mass
value of a meshless node is spread around the node with the kernel function.
Therefore the density of a meshless node is calculated after the mass allocation
step by taking the weighted average of the masses of the neighboring nodes
i m j w(rij )
(3)
u u X
uY
uZ x X .
T
u uT uuT
,
(4)
where the gradient of the continuous displacement vector field u is essentially
the derivatives of (ux , u y , uz ) with respect to ( X , Y , Z ) arranged in the Jacobian
format. For an isotropic linear-elastic material, the strain is mapped to the stress
by the C tensor that approximates the material properties and is composed of
two independent coefficients, Youngs Modulus and Poissons Ratio
C .
(5)
We need the partial derivatives of the displacement vector field in order to
compute the strain, stress, and the internal elastic forces applied to the meshless
particles. Moving Least Squares (MLS) approximation is used to compute this
gradient ( u ).
For a central meshless particle i and its neighbor j, the value of u x at the location of j can be approximated by the first order Taylor expansion as
ux j uxi
u X
.( X j Xi ) .
X i
(6)
The weighted sum of squared differences between the displacement vector and
its approximation obtained from the equation (6) gives the error measure of the
MLS approximation
e (ux j ux j )2 w(rij ) .
(7)
The expanded equation of the error measure for a particle i is therefore obtained by
e (uxi (
j
ux
u
u
) Xij ( x )Yij ( x ) Zij ux j ) 2 w(rij ) .
X
Y
Z
(8)
We want to minimize this error measure for some values of (ux)/X, (ux)/Y,
and (ux)/Z, therefore we set the derivative of the error measure e with respect
to the partial derivatives of the displacement vector to zero, resulting in three
equations for three unknowns
ux
T
ux j uxi ( X j Xi ) w rij . (9)
( X j Xi )( X j Xi ) w rij
X
j
j
i
The coefficient of the partial derivative on the left-hand side of the equation (9)
is the 3x3 matrix called the moment matrix (A). A can be inverted and premultiplied with both of the sides of the equation for computing the partial derivatives.
U i vi
1
( i i ) .
2
(10)
The elastic force per unit volume at a meshless nodes location is the negative
directional derivative of the above strain energy density with respect to this nodes
displacement. The forces applied to the particle i and its neighbors j are then
fi ux U i vi i ux i
i
f j ux U i vi i ux i
j
(11)
Eq. (11) is iterated over all particles and the total force applied to a particle is
the sum of all the forces with the same index as that particle. These force components are obtained by using the Green-Saint-Venant strain tensor, which measures
the linear and shear elongation. In case of a volume inverting displacement
though, this strain becomes zero. In order to introduce restoring body forces in
cases of volume inversion, Muller et al. [19] added another energy term to the system that penalizes deviations from a volume conserving transformation
Ui
1
kv ( J i 1) 2 .
2
(12)
In this energy term, J i is the Jacobian of the displacement vector field mapping and kv is the volume restoration constant. Although in our experiments we
have not worked on cases that cause the inversion of the volume, the effect of this
term and the volume restoration constant parameter in cases of large deformations
is yet to be examined.
For each of the meshless nodes, these force components are accumulated in the
force vectors and then passed to the SOFA time integration module. The modular
implementation of the SOFA library allows the use of different time integration
schemes without changing the actual implementation of the algorithm. In this
work, we used the 4th order Runge-Kutta integration scheme.
functionalities of the library by deriving new components from the existing ones.
The meshless elastic model described in this work has been implemented as a
force field component, which is easily interchangeable with the existing force field
components such as hexahedral finite elements or mass-spring networks. The required component interfaces such as initialization and force accumulation were
implemented according to the algorithm steps described in the previous section.
Due to the nature of the presented verification method, accurate contact handling is essential to obtain precise results. Unfortunately, because of the approximating nature of the meshless methods, imposing Dirichlet boundary conditions is
a challenging task on its own [24]. In this work, we have followed the approaches
that were originally adopted in the SOFA library. In SOFA, constraints are filter
like components, which cancel out the forces and displacements applied to their
associated particles. For example for fixed node boundary conditions, SOFA's
fixed constraint component is used to attach a meshless particle to a fixed point in
the current configuration.
SOFA has support for several contact handling methods such as penalty-based
and constraint-based methods. Among these, the constraint-based methods are
more appealing than the former class of methods because they use Lagrange multipliers to handle complex constraints and produce physically accurate results with
the additional computation cost. Lagrange multipliers with unilateral interaction
laws are used to handle complex constraints. The constraints depend on the relative positions of the interacting objects, which are the meshless particles and the
spherical rigid indenter in our case [12].
4 *
E Rd 3 2 ,
3
(13)
where F is the vertical force applied on the spherical load, R is the radius of the
1 1 v12 1 v22
.
E*
E1
E2
(14)
10
The Hertzian theory assumes 1) small strains within the elastic material, 2)
much smaller area of contact compared to the areas of the objects in contact, and
3) continuous and frictionless contact surfaces.
The Hertzian theory is an important stepping stone in the field of contact mechanics; therefore there have been numerous finite element analysis studies about
the subject that use both research and commercial finite element code [26-29].
Figure 1 Comparison of the FEBio FEM Code and the theoretical solution of
the Hertzian non-adhesive frictionless contact theory.
11
The contact mechanics experiment was setup as a SOFA scene. In order to assess the ground-truth performance of the contact handling in SOFA, the rectangular deformable block was represented with the hexahedral finite element model in
addition to the described meshless method. The validation of the hexahedral FEM
implementation of SOFA is studied by Marchal et al. [32]. For a given sphere radius, simulations were performed for varying force values (F) applied to the spherical load (Figure 2). With the applied force, the rigid sphere comes into contact
with the block and deforms it. The vertical velocity of the sphere is monitored and
the indentation of the material (d) is measured when the sphere comes to rest. This
F-d pair is compared to the theoretical solution.
(a)
(b)
Figure 2 (a) Initial setup of the indentation experiment for the SOFA FEM
model, (b) the close-up view of the indented deformable material. SOFA allows the user to track and monitor simulation values of indexed particles.
The meshless nodes are distributed uniformly inside a cubical volume with 2m
long edge length. The indentation experiment is repeated for several distribution
configurations, which play critical role especially for the MLS approximationbased collocation methods. The convergence rate in the L2 (vector) error norm of
the force-indentation pairs with respect to the theoretical values (Figure 3) is investigated. The effect of different distribution schemes on the accuracy, stability,
and performance of the meshless collocation methods is yet to be examined. In our
implementation, the number of neighboring nodes is limited to 16 for each of the
meshless nodes.
12
Figure 3 Error in the L2 norm with respect to the theoretical solution as functions of total number of the degrees of freedom for the (a) FEM model and
(b) meshless method.
We compared the SOFA FEM implementation and the meshless collocation
method with close accuracy (Figure 4). For our meshless collocation method with
nodal integration, we used an explicit time integration scheme with a time step of
0.001s without any stability problem. For the SOFA FEM implementation, we had
to use implicit integration with a time step of 0.01s. The calculations were performed within the SOFA application on a single Intel Core i5 CPU running at 2.67
GHz with 16 GB of RAM under Windows 7 operating system. The SOFA FEM
implementation took 195ms of calculations per time step, whereas the meshless
method consumed 20.11ms for calculations per time step. Therefore, the meshless
collocation implementation in SOFA (along with other SOFA related operations
such as collision detection) is roughly 25 times slower than the real-time operation, which is slightly better than the 30 times slower performance reported by the
Meshless TLED algorithm [20]. The calculation speed of the meshless collocation
algorithm is governed by the number of particles and the number of neighbors assigned to each particle.
13
14
References
[1]
[2]
[3]
[4]
[5]
[6]
[7]
[8]
[9]
[10]
[11]
[12]
[13]
[14]
[15]
[16]
[17]
15
[18]
[19]
[20]
[21]
[22]
[23]
[24]
[25]
[26]
[27]
[28]
[29]
[30]
[31]
[32]