OUTLINE
Models: General Overview
Mechanics and Continuum Mechanics
Mechanics of Solid Objects and Elasticity
Kinematics: displacements, deformations, strains
Kinetics: forces, pressures, stresses, tractions
Linear Elasticity: continuous formulation, FEM, solution
Hyperelasticity: towards non-linear models
Co-rotational approach: geometry-based compromise
MODELS
A model is an abstract structure that uses mathematical language to
describe the behaviour of a system.
typical examples of models:
– electrophysiological model: describes electrical properties of tissue
(e.g. electrophysiological model of heart)
– model of fluid dynamics: describes behaviour of liquid (e.g.
cardiovascular fluid mechanics (blood circulation)
– biomechanical model of an organ: describes elastic/plastic
behaviour of tissues (e.g. hyperelastic model of liver)
the mathematical language is usually based on differential equations
– the behaviour is “a change of state” (derivative)
MECHANICS
area of science dealing with physical bodies subject to force and/or
displacements
classical (Newtonian) vs. quantum mechanics :-)
– kinematics (geometry of motion): moving points/bodies without
considering the causes of motion
– (analytical) dynamics: relationship between motion of bodies and its causes
CONTINUUM MECHANICS
deals with the analysis of the kinematics and the mechanical behavior
of materials modeled as a continuous mass rather than as discrete
particles
continuum hypothesis: well defined properties in infinitely small
points (reference element of volume)
solid mechanics: study of continuous materials with defined rest shape
fluid mechanics: study of fluid materials (liquids, gases, plasmas)
e.g. CFD (computational fluid dynamics)
obeying common laws: conservation of mass, energy, [linear and
angular] momentum
SOLID MECHANICS
studies the behavior of solid materials, especially their motion and
deformation under the action of forces, temperature changes, phase
changes, and other external or internal agents.
elasticity: describes materials that return to their rest shape after
applied stresses are removed
viscoelasticity: elastic material with damping (hysteresis loop)
plasticity: describes materials that permanently deform after a
sufficient applied stress
thermoplasticity: coupling between mechanics and thermal properties.
ELASTICITY
ability of a body to resist a distorting influence or stress and to return to
its original size and shape when the stress is removed
basically, it defines mathematic relation between displacements and
applied forces
kinematics: relates displacement to strain (geometry)
kinetics: relates forces to stresses (e.g. equilibrium)
constitutive law: relation between the stress and strain (the material)
linear elasticity: keeping all relations linear (non-conservative!)
hypoelasticity: extension of linear elasticity
hyperelasticity: a family of models (materials), typically used for tissues
TOWARDS THE LINEAR
ELASTICITY
VECTOR AND TENSOR
FIELDS I
continuum mechanics: body as a continuum set of particles (3D points)
initial configuration X (X,Y,Z) vs. deformed configuration x (x,y,z)
displacement – vector function in 3D defined for in each particle (vector
field)





elasticity theory formulated using tensors
similarly as vector field, tensor field is a “tensorial” function defined
in each particle (i.e., over the continuous domain)
typical operators on fields: gradient, divergence, curl
u(x, y, z) = (ux(x, y, z), uy(x, y, z), uz(x, y, z))
x = X + u
VECTOR AND TENSOR
FIELDS II
Tensor notation:
–summation over repeated indices
–derivative using ‘,’ notation
aijbj ⌘
X
j
aijbj
fi,j ⌘
fi
xj
Vector-matrix notation:
–using bold symbols: A, σ (matrix), v (vector)
–derivatives written as operators: gradient: rf = ( @
@x , @
@y , @
@z )>
f
Example:
–divergence of a vector field
divu = (@ux
@x +
@uy
@y + @uz
@z ) = ( @
@x , @
@y , @
@z ) · (ux, uy, uz)>
= r · u = ui,i
u(x, y, z) = (ux(x, y, z), uy(x, y, z), uz(x, y, z))
Constitutive
equation
STATIC LINEAR
ELASTICITY
Kinematics Kinetics
KINEMATICS:
DEFORMATION
deformation field: vector field defined in each point 

deformation gradient: 2nd order tensor defined in each point

decomposition of deformation gradient to rotation and stretch
tensors
right Cauchy-Green deformation tensor (square of local change)

alternative: left Cauchy-Green deformation tensor
F = I + ru
C = F>
F = I + ru + ru>
+ ru>
ru
x = X + u(x, y, z)
F = RU = V R : R 1
= R>
B = FF>
= I + ru + ru>
+ ru>
ru
KINEMATICS: STRAIN
strain: a description of deformation in terms of relative displacement
of particles in the body that excludes rigid-body motions
different measures of strain: Green, Biot, Almansi, logarithmic strain
Green strain tensor:

linearization:








E =
1
2
(C I) =
1
2
(ru + ru>
+ ru>
ru)
geometric non-
linearity
" = e = 1
2 (ru + ru>
)
KINEMATICS: STRAIN
components of strain: diagonal + shear strains:



















Constitutive
equation
ELASTICITY-BASED
MODELING
Kinematics
Strain –
Displacement
Kinetics
" = 1
2 (ru + ru>
)
KINETICS: STRESS
stress: internal forces that neighboring particles of a continuous
material exert on each other
Cauchy (true) stress tensor: 2nd order tensor that completely define
stress at a point
relates a unit length vector and stress vector:
the components of stress vector (surface traction):
t = n
ti = dgi
dS
STRESS TENSOR
stress: internal forces that neighboring particles of a continuous
material exert on each other
Cauchy (true) stress tensor: 2nd order tensor that completely defines
stress at a point
conservation of linear momentum: in static equilibrium, it satisfies
equilibrium equation in each point (b being the body forces)

conservation of angular momentum: symmetry (6 components
instead of 9)
ij = ji
⌧xy = ⌧yx
⌧xz = ⌧zx
⌧yz = ⌧zy
div + b = 0 i.e., r · + b = 0 i.e., ij,j + bi = 0
Constitutive
equation
ELASTICITY-BASED
MODELING
Kinematics
Strain –
Displacement
Kinetics
Stress in static
equilibrium
" = 1
2 (ru + ru>
) r · + b = 0
t = n
CONSTITUTIVE
EQUATION
Cauchy elastic material: stress is a function of strain
linear elasticity: stress is a linear function of strain
Hooke law: the relation between stress (2nd order tensor) and strain
(2nd order tensor) is a 4th order tensor

in general, C has 81 components: however, symmetry of strain and
stress reduces the number of components to 21
for isotropic and homogeneous material, number of parameters is
reduced to two Lamé coefficients:
= Itr(") + 2µ"
ij = Cijkl"kl i.e., = C : "
MATERIAL PARAMETERS
in tensorial notation (with Einstein summation convention):

Lamé coefficients: the second is sometimes called shear modulus (G) 



where
E is the Young’s modulus [Pa]: stiffness of the material
nu is the Poisson’s ratio: incompressibility of the material <0,0.5
= Itr(") + 2µ"
ij = ij"kk + 2µ"ij = ijuk,k + µ(ui,j + uj,i)
C100 = µ C200 =
⌅ + 2µ
8
C010 =
µ
3
(2.2
ere µ and ⌅ are Lamé material constants. The energy function W is then given as follow
W =
⌅
2
ii jj + µ ij ji. (2.2
ng the relation 2.20, the stress/strain relation is linear:
ij = ⌅ mm⇥ij + 2µ ij (2.2
refore, the only non-linearity in the StVenant material is given by the displacement/str
tionship shown by the Eq. 2.7. The two Lamé coefficients from the Eq. 2.23 are define
⌅ =
E⌃
(1 + ⌃)(1 2⌃)
µ =
E
2 + 2⌃
(2.2
ere E is Young modulus and nu ⇥ (0, 1
2) is the Poisson ratio which determines th
ompressibility of the material. The values of both E and ⌃ has been determined exper
ntally for various types of materials and they can be found in literature.
Second, in the case of Mooney-Rivlin material, Cijk = 0 except for C100 and C010. Th
ed energy function W is determined as
⇥
Constitutive
equation
ELASTICITY-BASED
MODELING
Kinematics
Strain –
Displacement
Kinetics
Stress, static
equilibrium
= Itr(") + 2µ"
Stress-strain relation
" = 1
2 (ru + ru>
) r · + b = 0
t = n
Navier-Cauchy equation (see the proof performed by components
on LinearElasticity@Wikipedia):
PUTTING IT ALL
TOGETHER
" = 1
2 (ru + ru>
) = Itr(") + 2µ"
tensor notation:
per component: K 2 {x, y, z}
r · + b = 0
t = n
( + µ) @
@K
⇣
@ux
@x +
@uy
@y + @uz
@z
⌘
+ µ
⇣
@2
uK
@x2 + @2
uK
@2y2 + @2
uK
@z2
⌘
+ bK = 0
( + µ)r(r · u) + µr2
u + b = 0
( + µ)uj,ij + µui,jj + bi = 0
the body given by a continuous domain with boundary
Navier-Cauchy equation holds for every point of the domain 

(fi being body forces per unit volume)



essential boundary conditions has to be defined on a part of the
boundary (to choose the particular solution of N.-C. PDE 

natural boundary conditions can be defined on a part of the
boundary (i.e., tractions T along normal n in point p)

THE PROBLEM TO SOLVE
˜⌦
Tp
i = ijnp
j for p 2 ˜N where ˜N ⇢ ˜ and ˜ = @ ˜⌦
up
i = ¯up
i for p 2 ˜E where ˜E ⇢ ˜ and ˜ = @ ˜⌦
˜ = @ ˜⌦
( + µ)uj,ij + µui,jj + bi = 0
the only feasible way – discretization: approximate the original
continuous quantities by discrete (piecewise) functions:
central role of the interpolation (basis, shape, blending) functions
required properties:
local support:

the function is non-zero

only inside the element
bound to a node n:

CONTINUOUS VS.
DISCRETE SOLUTION II
CHAPTER 2. THEORETICAL BACKGROUND 13
x
0.8
-0.4
0.5 1
0.4
-0.5 0-1
0
(a)
-1
-0.5
y0
0.5
1
-1
-0.5
0
0.5 x0
1
0.2
0.4
0.6
0.8
1
(b)
-1
-0.5
y
0
0.5
1-1
-0.5
x
0
0.5
1
0
0.1
0.2
0.3
0.4
0.5
(c)
Figure 2.2: Illustration of shape functions over finite elements. Both linear and quadratic functions
over one-dimensional domain are shown by (a). Linear shape function over 2D rectangular element
is depicted by (b) and the same 2D element is used with quadratic shape function at (c).
u(x) ⇡
X
n
Un'n(x)
@u(x)
@x
⇡
X
n
Un
@'n(x)
@x
'n(xm) = nm
FINITE ELEMENT
METHOD
First appeared in 40s and 50s (civil engineering, aeronautics).
1. Weak formulation of the continuous differential problem
– integration over domain and multiplication by test functions
2. Discretization 

– discretization of the domain by the elements
– discretization of the variable and the operator
– integration over element volume (quadratures)
3. Global assembling of the algebraic system of equations

– imposing the compatibility between the elements
4. Imposition of the essential boundary conditions
5. Numerical solution of the algebraic system
EXAMPLE: STATIC LINEAR
ELASTICITY (SLE)
ij,j + bi = 0 eij =
1
2
(ui,j + uj,i) ⇤ij = ⇥ekk ij + 2µeij
Given relations (in tensor notation)
Newton’s law (kinetics) linearized strain (kinetics) linear material (constitutional law)
ij,j + bi = 0 eij =
1
2
(ui,j + uj,i) ⇤ij = ⇥ekk ij + 2µeij
Given relations (in tensor notation)
Newton’s law (kinetics) linearized strain (kinetics) linear material (constitutional law)
Weak form of the Newton’s equation (Lax-Milgram lemma)
–integration over the volume 

–multiplication by a test functions wi
The integral over volume allows to distribute the derivatives
–application of chain rule

–divergence theorem
– no derivative of the stress tensor 

– the only derivative applied to the test function on the left side
– ti: tractions defined over the surface (natural boundary conditions)@⌦
Z
⌦
( ij,j + bi)wid = 0
Z
⌦
ijwi,jd⇥ =
Z
⌦
biwid⇥ +
Z
⌦
tiwid
SLE: DISCRETIZATION
AND GALERKIN METHOD
Galerkin method: use the same interpolation functions to discretize the test
functions w and the solution u over an element e:
Z
⌦
ijwi,jd⇥ =
Z
⌦
biwid⇥ +
Z
⌦
tiwidThe actual weak form:
Domain discretization by elements e:

– element e given by N nodes

– each element “equipped” with interpolation functions 

– index n: node of the element (therefore N interpolation functions per element)
where: ⇤ij = ⇥ekk ij + 2µeij eij =
1
2
(ui,j + uj,i)
en
(x, y, z)
wi = en
Wen
i ui = en
Uen
i
Example of derivative:
(note: no summation over e!)
˜⌦ ⇡ ⌦ =
U
e ⌦e
wi,j = 'en
,j Wen
i
SLE: GALERKIN METHOD
II
Discretized week form:
where: ⇤ij = ⇥ekk ij + 2µeij eij =
1
2
( en
,j Uen
i + en
,i Uen
j )
Galerkin method: the equations hold for any virtual displacement Wi:
X
e
Z
⌦e
ij⇥en
,j d⇥ Wen
i =
X
e
✓Z
⌦e
bi⇥en
d⇥ +
Z
⌦e
ti⇥en
d
◆
Wen
i
X
e
Z
⌦e
ij⇥en
,j Wen
i d⇥ =
X
e
Z
⌦e
bi⇥en
Wen
i d⇥ +
Z
⌦e
ti⇥en
Wen
i d
For each element e, we have the local equation:
Z
⌦e
ij⇥en
,j d⇥ =
Z
⌦e
bi⇥en
d⇥ +
Z
⌦e
ti⇥en
d
where: ⇤ij = ⇥⌅ne
,k Une
k ij + µ(⌅en
,j Uen
i + ⌅en
,i Uen
j )
SLE: THE ELEMENT
EQUATION
Z
⌦e
ij⇥en
,j d⇥ =
Z
⌦e
bi⇥en
d⇥ +
Z
⌦e
ti⇥en
d
where: ⇤ij = ⇥⌅ne
,k Une
k ij + µ(⌅en
,j Uen
i + ⌅en
,i Uen
j )
Right-hand side:
– we consider tractions to be zero and
– body forces to be constant w.r.t. space
Left-hand side:
– clearly linear in U being the unknown displacements in nodes n=1...N
Z
⌦e
⇥⇤ne
,k Une
k ij + µ(⇤en
,j Uen
i + ⇤en
,i Uen
j )⇤en
,j d
– since linear, the left-hand side can be re-organized to Ken
ij Uen
j
bi
Z
⌦e
'ne
d⌦
VOIGT NOTATION
Left-hand side:
– the tensor notation has been useful to derive the final form
– for implementation purposes, Voigt notation is usually employed where 3x3
symmetric 1-order tensor is stored as 6x1 vector:
Z
⌦e
ij⇥n
,jd⌦
⇤ij = ⇥ekk ij + 2µeij
eij =
1
2
( n
,jUn
i + n
,iUn
j )
with
SLE: STRESS-STRAIN
MATRIX D
Applying the Voigt notation to the stress–strain relation

results in following matrix equation (derivation is straightforward:

The matrix in the middle is 6x6 stress-strain matrix (denoted further as D).

Before encoding the rest into matrices we have to choose the interpolation functions!
Z
⌦e
ij⇥n
,jd⌦
Note that only derivatives of interpolation functions appear in the formulation.
0
B
B
B
B
B
B
@
⇥11
⇥22
⇥33
⇥12
⇥13
⇥23
1
C
C
C
C
C
C
A
=
0
B
B
B
B
B
B
@
+ 2µ 0 0 0
+ 2µ 0 0 0
+ 2µ 0 0 0
0 0 0 µ 0 0
0 0 0 0 µ 0
0 0 0 0 0 µ
1
C
C
C
C
C
C
A
0
B
B
B
B
B
B
@
e11
e22
e33
2e12
2e13
2e23
1
C
C
C
C
C
C
A
⇤ij = ⇥ekk ij + 2µeij
eij =
1
2
( en
,j Uen
i + en
,i Uen
j )
P1: TETRAHEDRAL
LINEAR ELEMENT
– tetrahedral: simplex in 3D having
four nodes
– linear since we choose linear
interpolation functions:
(a general linear function in 3D)
– how to find the coefficients a,b,c,d? Recall the basic property of an
interpolation function:
(the value of an interpolation function associated to a node i is 1 when
evaluated in that node [xi, yi, zi] and zero in any other node [xj,yj,zj])
⇥i
(xj, yj, zj) = ij i, j 2 1, . . . N
(x, y, z) = a + b(x) + c(y) + d(z)
SLE&P1: COMPUTING THE
SHAPE FUNCTIONS
Linear P1 (Lagrangian) tetrahedral element
– putting the condition into a matrix form gives:
– denoting V the matrix on the left (nodal matrix), 4 instances of coefficients
corresponding to 4 interpolation functions (associated to each node) can be
computed as columns of the V–1 (recall the requirements for mesh quality!)
– recall also that only derivatives of interpolation functions are present in the
formulation (so only coefficients b,c,d) will be used
0
B
B
@
1 x1 y1 z1
1 x2 y2 z2
1 x3 y3 z3
1 x4 y4 z4
1
C
C
A
0
B
B
@
a
b
c
d
1
C
C
A =
0
B
B
@
1 0 0 0
0 1 0 0
0 0 1 0
0 0 0 1
1
C
C
A
SLE&P1: STRAINDISPLACEMENT
MATRIX B
Using the Voigt notation and assuming the linear P1 tetrahedra used for
discretization, the left-hand side
can be rewritten in matrix form as:
Z
⌦e
ij⇥n
,jd⌦
De =
0
B
B
B
B
B
B
@
+ 2µ 0 0 0
+ 2µ 0 0 0
+ 2µ 0 0 0
0 0 0 µ 0 0
0 0 0 0 µ 0
0 0 0 0 0 µ
1
C
C
C
C
C
C
A
Z
⌦e
B>
e DeBed⌦
Be =
0
B
B
B
B
B
B
@
b1 0 0 b2 0 0 b3 0 0 b4 0 0
0 c1 0 0 c2 0 0 c3 0 0 c4 0
0 0 d1 0 0 d2 0 0 d3 0 0 d4
c1 b1 0 c2 b2 0 c3 b3 0 c4 b4 0
d1 0 b1 d2 0 b2 d3 0 b3 d4 0 b4
0 d1 c1 0 d2 c2 0 d3 c3 0 d4 c4
1
C
C
C
C
C
C
A
⇤ij = ⇥ekk ij + 2µeij eij =
1
2
( en
,j Uen
i + en
,i Uen
j )
SLE&P1: LOCAL STIFFNESS
MATRIX
What about the integration?
– recall that only derivatives of shape functions appear in the formulation
– since interpolation functions are linear, only coefficients b,c,d appear in the
matrices
– therefore, the integrand is constant (does not depend on x,y,z)
– integration of a constant over a tetrahedron is computed by multiplication of
the constant by the volume of the tetrahedron
– the volume of a tetrahedron is given by determinant of nodal matrix:

– the final form is therefore:
– the local matrices Ke are assembled into a global matrix K
– the contribution from different elements to the same node are added
(globalization matrix)
Ke =
Z
⌦e
B>
e DeBed =
|Ve|
6
B>
e DeBe
Ve =
|Ve|
6
ASSEMBLING THE
GLOBAL SYSTEM
the procedure now gives 12x12 matrix (4x4 block matrix where each
block (i,j) corresponds to stiffness relation between nodes n and m
(n,m=1…4)
global assembly:
mapping for each node from local to global indices: (e,n) -> n
the block (n,m) from matrix associated to element e is added to
the global block at position (n,m) in the global matrix
usually is done directly during the computation of local matrix
the global matrix is a 3Nx3N block matrix where N is the total
number of DOFs (and 3N is thus the number of degrees of
freedom)
BOUNDARY CONDITIONS
choosing a particular solution (otherwise K singular)
several options to impose a Dirichlet boundary condition ui=V
elimination (projection):

— left side: K(i,k) = K(k,i) = 0 for all k ≠ i, K(i,i) = 1

— right side: f(i) = V (“pseudo-loads”)

— not very flexible and difficult to parallelize
penalization: adding a penalization term to impose the boundary
condition (reduces the “quality” of matrix in terms of the
condition number)
Lagrange multipliers: changes the properties of the matrix
(larger, possibly indefinite)
THE GLOBAL STIFFNESS
MATRIX
linear relation between forces (f) and displacements (u):
encoding relations between nodes
highly sparse (<3% of non-zero)
non-zero blocks only for combinations

of nodes connected by a mesh edge
suitable representation [i j Kij]
efficient matrix–vector multiplication
regular after the imposition of boundary conditions
symmetric, positive-definite, sparsity pattern depends on node
numbering (can be improved e.g. by Metis)
Ku = fCHAPTER 2. THEORETICAL BAC
(a)
Figure 2.3: Illustration of the sparsity
straints applied (b) boundary condition
and boundary conditions applied via L
PRACTICAL MATRIX
MANIPULATION
sparse matrices generated from the FE formulation
only a small fraction of entries non-zero (<3%)
system of N nodes in 3D results in size of (3N)2
practical example: 10000 nodes in double (4B): 3.4GB
but 3.3GB are zeros...
common format: i j Aij (137MB, 2 x int + 1 x double)
row vs. column compressed
sometimes storing both representations can be practical
SYSTEMS OF LINEAR
EQUATIONS
scalar case: ax = b x =
b
a
vectorial case: Ax = b x = A 1
b
properties of A (considered being a square matrix)
regular matrix: inverse A-1 exists
symmetric: equals to transpose, AT = A
positive-definite: zTAz is positive for a vector z (eigenvalues)
orthogonal matrix: AT = A-1 (representation of rotations)
DIRECT SOLUTION OF
LINEAR SYSTEM
solution x = A 1
b
direct solutions: the inverse A–1 computed explicitly as factorization
for cases when you need to recompute Ax=b’ for another b’
2 phases: decomposition (factorisation), solution (back-substitution)
Cholesky decomposition: A = LLT (L lower triangular matrix): symmetric
positive-definite matrices, most optimal (num. of operation)
LDL decomposition: A = LDLT (D diagonal), works for some indefinite
matrices where Cholesky fails
LU decomposition: (U upper triangular matrix), general case, modified
Gaussian elimination (Doolittle, Crout algorihms, pivoting)
ITERATIVE SOLUTION OF
LINEAR SYSTEM
solution
will depend on properties of A
x = A 1
b
iterative solutions: the inverse A–1 is not assembled explicitly
start with an estimation x(0) and iterate until |Ax(i)–b| < e (stopping
criterium usually more complicated, absolute vs. relative residual )
conjugate gradients (CG): for symmetric, positive-definite matrices
(see Shewchuk: Conjugate gradients without agonizing pain)
bi-conjugate gradient (BiCG): generalization for non-symmetric
generalized minimal residual (GMRES): any regular matrix
preconditioned versions: approximation of A–1
ISSUES WITH LINEAR
ELASTICITY
after imposition of the boundary conditions, 

the system can be solved
iterative: even the matrix K does not have to be assembled
direct: the both K and K–1 are assembled and stored explicitly,
so u can be updated for any new f
CHAPTER 2. THEORETICAL BACKGROUND 8
(a) non-linear (b) non-linear (c) linear (d) linear
Figure 2.1: Illustration of the difference between the models employing non-linear strain tensor
(figures (a) and (b)) and linear strain tensor L
(figures (c) and (d)). The figure (d) is scaled by factor
CHAPTER 2. THEORETICAL BACKGROUND 8
(a) non-linear (b) non-linear (c) linear (d) linear
Figure 2.1: Illustration of the difference between the models employing non-linear strain tensor
(figures (a) and (b)) and linear strain tensor L
(figures (c) and (d)). The figure (d) is scaled by factor
CHAPTER 2. THEORETICAL BACKGROUND 8
(a) non-linear (b) non-linear (c) linear (d) linear
Figure 2.1: Illustration of the difference between the models employing non-linear strain tensor
(figures (a) and (b)) and linear strain tensor L
(figures (c) and (d)). The figure (d) is scaled by factor
EORETICAL BACKGROUND
(b) non-linear (c) linear (d) linearlinearized Green strain does not work for large deformations
Ku = f
TOWARDS NONLINEAR: COROTATIONAL
FORMULATION
an extremely successful approach in soft-tissue modeling allowing
for large displacements (but supposing small strains)
C.Felippa: A systematic approach to the element-independent corotational dynamics of
finite elements, 2000
uses the linear-elasticity 

but co-rotational strain
the simulation is performed 

in small steps and in each step:
the actual deformation of every element e is decomposed into rigid and
deformable components w.r.t. the initial configuration
the rigid component is given by a rotation Re of the component
the local stiffness matrix Ke is updated as R>
e KeRe
Matthieu Nesme, Yohan Payan & Franç
matrix: J = e0
1 e0
2 e0
3
−1
[e1 e2 e3], with respect to the
tial state, where the e0
i are the initial edge vectors and th
are the current ones. Matrix J is then decomposed in o
to extract separately a rigid rotation R applied to the elem
and a deformation E as shown fig. 1. This decompositio
not unique and several approaches can be considered.
in the initial local frame
at rest form
deformed and replaced
R
E
J=RE
deformed
and displaced
Figure 1: An initial tetrahedron is deformed by the trans
mation J composed both a rigid motion R and the defor
tions are contained in E.
Polar decomposition Etzmuß et al [EKS03], followed
other authors [HS04, MG04], presented a method ba
on the polar decomposition using eigenvalues and eig
vectors. The polar decomposition of a square matrix c
putes the nearest orthogonal frame to the given colu
CO-ROTATIONAL
FORMULATION II
the matrix K is not constant anymore ( K => K(u) )
the rotational matrices Re(u) depend on the actual u
in each step, Newton-Raphson method should be performed, actually,
works quite stably even if only one iteration is performed
the decomposition can be performed by various methods
choosing the basis
polar(1), QR(2), SVD
although the large deformations are simulated realistically, only
small strains are handled correctly
more information about the implementation in SOFA:
M.Nesme et al.: Efficient, physically plausible finite elements, 2005
sme, Yohan Payan & François Faure / Efficient, physically plausible finite elements
2 e3], with respect to the iniitial
edge vectors and the ei
s then decomposed in order
ion R applied to the element
g. 1. This decomposition is
es can be considered.
deformed and replaced
R
b
c
z
a
y
x
d
(1)
c
b
d z
x
y
a
(2)
Figure 2: the local frames. (1) Polar decomposition: singl
frame reflecting best the matter, nearest to the edges. (2) QR
decomposition: the first axis is the first edge ab, the second
axis is orthogonal to the first on plane (ab,ac), and the las
axis is obtained by construction of an orthonormal frame.
2.2. Newton’s law
Newton’s law on linear acceleration relates the acceleration
Nesme, Yohan Payan & François Faure / Efficient, physically plausible finite elements
e2 e3], with respect to the iniinitial
edge vectors and the ei
is then decomposed in order
ation R applied to the element
fig. 1. This decomposition is
hes can be considered.
deformed and replaced
R
b
c
z
a
y
x
d
(1)
c
b
d z
x
y
a
(2)
Figure 2: the local frames. (1) Polar decomposition: singl
frame reflecting best the matter, nearest to the edges. (2) QR
decomposition: the first axis is the first edge ab, the secon
axis is orthogonal to the first on plane (ab,ac), and the las
axis is obtained by construction of an orthonormal frame.
2.2. Newton’s law
Newton’s law on linear acceleration relates the acceleratio