Department of Science and Technology
Institutionen för teknik och naturvetenskap
Linköping University Linköpings Universitet
SE-601 74 Norrköping, Sweden
601 74 Norrköping
LiU-ITN-TEK-A--09/008--SE
Implementation and
Modification of a Cloth
Simulator for Deformable
Shells
Linus Hilding
Implementation and
Modification of a Cloth
Simulator for Deformable
Shells
Examensarbete utfört i medieteknik
vid Tekniska Högskolan vid
Linköpings universitet
Linus Hilding
Handledare Luigi Tramontana
Examinator Matt Cooper
Upphovsrätt
Detta dokument hålls tillgängligt på Internet – eller dess framtida ersättare –
under en längre tid från publiceringsdatum under förutsättning att inga
extra-ordinära omständigheter uppstår.
Tillgång till dokumentet innebär tillstånd för var och en att läsa, ladda ner,
skriva ut enstaka kopior för enskilt bruk och att använda det oförändrat för
ickekommersiell forskning och för undervisning. Överföring av upphovsrätten
vid en senare tidpunkt kan inte upphäva detta tillstånd. All annan användning av
dokumentet kräver upphovsmannens medgivande. För att garantera äktheten,
säkerheten och tillgängligheten finns det lösningar av teknisk och administrativ
art.
Upphovsmannens ideella rätt innefattar rätt att bli nämnd som upphovsman i
den omfattning som god sed kräver vid användning av dokumentet på ovan
beskrivna sätt samt skydd mot att dokumentet ändras eller presenteras i sådan
form eller i sådant sammanhang som är kränkande för upphovsmannens litterära
eller konstnärliga anseende eller egenart.
För ytterligare information om Linköping University Electronic Press se
förlagets hemsida
http://www.ep.liu.se/
Copyright
The publishers will keep this document online on the Internet - or its possible
replacement - for a considerable time from the date of publication barring
exceptional circumstances.
The online availability of the document implies a permanent permission for
anyone to read, to download, to print out single copies for your own use and to
use it unchanged for any non-commercial research and educational purpose.
Subsequent transfers of copyright cannot revoke this permission. All other uses
of the document are conditional on the consent of the copyright owner. The
publisher has taken technical and administrative measures to assure authenticity,
security and accessibility.
According to intellectual property law the author has the right to be
mentioned when his/her work is accessed as described above and to be protected
against infringement.
For additional information about the Linköping University Electronic Press
and its procedures for publication and for assurance of document integrity,
please refer to its WWW home page:
http://www.ep.liu.se/
the bending energy model. Recent articles have highlighted the close connection between the physical behaviour of cloth and that of thin plates and shells. By altering the bending model, the application is able to create animations of surfaces with higher bending stiff-ness than cloth simulators in general, at interactive frame rates. It also makes it possible to handle thin shells - surfaces that are curved in the undeformed state. Another feature looked for when deciding on a model, is the ability to simulate surfaces than do not stretch or shear significantly. This commonly means a very high stiffness coefficient for the forces. This is handled in the implementation by using a semi-implicit time stepping scheme. The simulation of plastic materials is briefly examined as well, by implementing a simple model for inelastic deformation. The work was carried out at the company Craft Animations in Gothenburg who, in their development of tools for physically based animation, may have an interest in the ability to simulate deformable thin surfaces.
Contents
Abstract i List of Figures iv 1 Introduction 1 1.1 Background . . . 1 1.2 Project description . . . 1 1.3 Thesis outline . . . 2 1.4 Limitations . . . 2 2 Theoretical background 3 2.1 Overview . . . 33 Implementation of a semi-implicit cloth simulator 5 3.1 Overview . . . 5
3.2 Geometric representation . . . 5
3.3 Backward Euler scheme . . . 6
3.4 Forces and their gradients . . . 7
3.4.1 Potential energy and condition functions . . . 7
3.4.2 The bend condition function . . . 8
3.4.3 The stretch condition function . . . 8
3.4.4 The shear condition function . . . 8
3.4.5 Differentiating the condition functions . . . 9
3.4.5.1 Finite differences . . . 9
3.4.5.2 Analytic expressions . . . 10
3.5 Dissipative forces . . . 11
3.6 Solving the system of linear equations . . . 12
3.6.1 Constraints . . . 12
3.6.2 Ascher and Boxerman’s modified pre-conditioned conjugate gradient method . . . 14
3.6.3 Exploiting the matrix sparsity . . . 15
3.7 Collision handling . . . 16
3.8 The simulation loop . . . 16
3.9 Conclusion . . . 16
4 Modification of the bending model 18 4.1 Quadratic Curvature Energies . . . 18
4.1.1 Implementation . . . 19
4.1.1.1 Results . . . 20
4.2 Cubic shells . . . 21
4.2.1 Implementation . . . 22
4.2.1.1 Changing in-plane energy model . . . 22 ii
List of Figures
3.1 Snapshot from a simulation in which a piece of cloth is completely constrained along one border edge. Gravity and wind are present. . . 15 3.2 Snapshots from a simulation. A piece of cloth, constrained along one edge,
falling over a solid immobile sphere. The mesh consists of 200 triangles and 121 vertices. The duration of the sequence is 55 seconds, with a frame rate of 8-10 fps on a 2.2 GHz Intel Core 2 Duo CPU. . . 17 4.1 Snapshots from a simulation, spanning over 10 seconds. A thin patch of a
semi-stiff material is falling down over a solid sphere. The surface will seek its rest state after impact. The mesh consists of 200 triangles and 121 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 26-28 fps. . 20 4.2 Visualization of the forces resulting from the cubic bending model, acting on
the vertices in a two-triangle hinge, having a rest angle of π/9. The angles have been deformed to π/4 and 3π/4 radians respectively, for the left and right image. The non-flat rest state is shown by the dashed green lines, and the flat state is shown by the dashed blue lines. The solid blue lines represent the forces acting to achieve a flat state, whereas the solid red lines depict the counter-acting non-flat forces. The resultant forces, shown by solid green lines, will act to achieve the non-flat rest state. . . 23 4.3 Snapshots from a simulation, spanning over 1.5 seconds. A thin shell object,
shaped like a bowl, is falling to a solid half-sphere, fixed to the ground. The surface will seek its curved rest state after impact. The mesh consists of 130 triangles and 77 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 62-64 fps. . . 25 4.4 Snapshots from a simulation with plasticity, spanning over 0.2 seconds. A
thin shell object, shaped like a bowl, is falling to a solid half-sphere, fixed to the ground. Notice how the dent is not fully restored after impact. The mesh consists of 130 triangles and 77 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 62-64 fps. . . 26
Introduction
1.1
Background
In computer graphics, the need for animations of moving cloth arises in areas such as film production, computer games and virtual reality applications. The physically based simula-tion of cloth and fabrics has been a research area for several decades. Another area with several similarities is that of thin plates and shells - thin walled structures, or 2D surfaces in 3D. Even though the latter has seen applications mostly in the field of structural engi-neering, where the computation models have been complex and computationally expensive, recent articles have explored methods that are simpler and better suited for use in computer graphics, where demands on accuracy are relaxed. The close connection between the phys-ical behaviour of cloth and that of thin plates and shells makes it possible to use similar techniques to simulate such diverse phenomenon as clothing, leaves, car bodies, and hats. The idea for the project came from the company Craft Animations being interested in investigating how a cloth simulator can be extended to be able to animate thin shells. Craft Animations is developing plug-ins for 3D Studio Max and Maya, providing the possibility to incorporate physically based simulation into the creation of animated scenes. Their main area of focus is the simulation of vehicles of various sorts. The possibility to handle deformable surfaces could be of interest for this type of application.
1.2
Project description
The underlying idea behind this thesis is to investigate how a cloth simulator framework can be extended with a novel bending model, to function as a simulator for thin plates and thin shells. With these features in focus, the demands on the cloth simulator regarding realistic reproduction of folds and wrinkles are relaxed. (This is otherwise a critical part of a realistic simulation of cloth, e.g. when simulating virtual garments.)
The thesis project consists mainly of two parts. The first is the implementation of a cloth simulator that is suitable for modifications with recent results from the field of discrete differential geometry. The second part is the implementation of these modifications.
Chapter 1. Introduction 2
1.3
Thesis outline
The next chapter of the thesis will be looking at some previous work within the field of cloth simulation. The aim of this section is to find a suitable technique for experimentation with the bend energy formulations that are often used for simulating the internal dynamics of cloth. Chapter 3 will describe the theory and implementation of a chosen cloth simulator framework. Chapter 4 will describe how some more recent formulations for the bending energy of a thin surface can be incorporated into the existing solution, and in what way it affects the behaviour and performance of the simulation. It will also describe how plasticity is added to the simulation.
1.4
Limitations
The cloth simulator framework chosen as theoretical foundation for the first stage of the implementation involved taking second order derivatives of energy functions [1]. These calculations are not explicitly derived in the article that was being followed, and experiments with this occupied more time than what was planned before starting the work. This resulted in lesser time being left for implementing a satisfying collision handling code, which was initially planned. Continuous collision detection is implemented by following the approach presented in [2] together with [3]. However, there was not enough time to get the response code to function properly. Therefore it was decided to omit this part of the simulator and leave it as an area of future development.
Theoretical background
2.1
Overview
The true nature of the internal dynamics of a fabric, with interacting woven fibres on the low level is often far too complex to model. Instead, the common approach is to treat the cloth as a 2-dimensional surface that resists deformations out of the plane and, to a far greater extent, in the plane. Some of the earlier experiments were carried out trying to predict the final drape of a cloth. In the field of Finite Elements, cloth has been subject to nonlinear continuum mechanics formulations that have managed to capture the behaviour quite accurately. However, these methods are often very computationally expensive, and they are considered challenging to implement.
To ease the implementation, and speed up the simulations, several simplified approaches have been tried within the graphics community. A common technique is the so called mass-spring model, in which a grid of particles, connected together by damped springs, represents a planar cloth surface. The particles are assigned a point mass, and when each particle is connected to its direct neighbours, the springs will act to prevent the particles from moving apart in the plane of the grid. If additional springs, connecting every particle to its second neighbour are introduced, they would work to prevent the cloth particles from moving out of the plane. This type of deformation is often referred to as bending within the cloth. Although the mass-spring model is simple to implement, it is hard to define a good model for non-uniform grids. And due to the non-physical derivation of the model, it is difficult to find appropriate stiffness coefficients for the spring forces. This drawback can be decisive when choosing an appropriate model, since the formulation for bend energy in a cloth simulation has significant impact on the behaviour of features such as wrinkles and folds. In addition, this technique can only be used to simulate surfaces with a flat rest state, i.e. thin plates, since a correct energy minimum of curved surfaces, i.e. thin shells, cannot be captured by a model of springs between second neighbours. The model is insensitive to the sign of the angle between two faces, and so called ”popping” may occur when large deformations suddenly change the sign of the angle at the rest state.
Over the past decade, there have been a number of articles trying to formulate various models for the bend energy. While the articles differ regarding the numerical solution techniques, many of them are similar in the sense that they formulate the bend energy as a scalar-valued function over the cloth surface, describing its potential energy. The forces that arise from this setting are the negative gradient of the potential energy, i.e. they act to achieve a state of equilibrium. Rather than using a global energy function
Chapter 2. Theoretical background 4 for the whole domain, it is often easier to use a locally defined function, discretized to a subset of the underlying representation of a particle grid. These subsets can then be evaluated separately, and then be assembled together to yield a representation of the entire system. In this sense the methods resemble that of the Finite Element Method (FEM), but are generally not as complex to implement. One main difficulty arises when using an implicit time discretization scheme, which necessitates taking second order derivatives of the corresponding energy function with respect to vertex degrees of freedom (DOFs), as will be discussed later on in the thesis.
There are several possible approaches to choose from when writing software to calculate derivatives of a function: Deriving an expression by hand and implement it as an expression, or using software with symbolic differentiation like Mathematica (which is often expensive). Although an analytic expression will give the exact derivative, the work of deriving the expression can indeed be tedious and error prone. One can also use the technique called Automatic Differentiation (sometimes referred to as Algorithmic Differentiation), which is a method for automatically taking derivatives of an algorithm by repeated use of the chain rule, yielding an exact result as well (down to machine precision). The derivatives can also be estimated using finite differences.
When using an explicit time integration scheme, this becomes much simpler. However, as real world fabrics tend to have a much stronger resistance to stretch than to bend or shear, the stretch forces need to be very strong. This results in a stiff equation system that generally requires small time steps. To improve performance, cloth simulators using an explicit integration method often release the demands of high stiffness coefficients for the stretch forces. Unfortunately, this results in a loose and bouncy, rubbery looking cloth. For the application of surfaces other than cloth, this is especially unacceptable. Some authors have however suggested an approach where bending forces are treated explicitly, and the stretch and shear forces are treated implicitly (see e.g. [4]).
In 1998, Baraff and Witkin published the paper Large steps in cloth simulation in which they present a full framework for a cloth simulator [1]. The major contribution of the paper is that it introduces a semi-implicit numerical integration method, which allows for using larger time steps than explicit methods do in general. The risk of instability still exists, but can easily be avoided by using a variable time step, e.g. letting it depend on the strain rate in the cloth. This makes it easier to maintain the desired inextensibility without degrading performance. Several authors have since contributed to improve the performance of their time discretization scheme, for example by addressing the unwanted damping that it suffers from, but it is nevertheless considered a seminal paper in the field. The article derives the internal forces of the cloth using a continuum mechanics formulation, discretized as per-triangle and per-hinge (two adjacent triangles, bending around their common edge) quantities for in-plane and out-of-plane motions respectively.
Due to the reported stability for large time steps, even with high membrane stiffness, the article of Baraff & Witkin is chosen as a good candidate for extension with the new bending models mentioned above. The next chapter will describe the implementation of a simulator based mainly upon this framework.
Implementation of a
semi-implicit cloth simulator
3.1
Overview
The following sections describe the implementation of the cloth simulator. At large, it is following the approach presented by Baraff and Witkin in the paper Large steps in cloth simulation [1]. In this article, they describe a full cloth simulation framework, complete with a backward Euler time integration technique and internal cloth dynamics models. They also provide a way to incorporate geometric constraints directly into the iterative solving of the resulting system of equations. The implementation of this thesis project uses their formulations for the internal cloth dynamics, as well as their time integration scheme. For solving the system of equations it uses the modified preconditioned conjugate gradient method presented by Ascher and Boxerman in [5].
Implementation of the cloth simulator is being done with the aim of keeping the physical simulation calculations and the geometric structure as separated as possible. This makes it easier to change one part of the structure without affecting the other.
3.2
Geometric representation
The cloth model adopted here is based on a system of particles connected together as a flat triangular mesh. A polygon mesh is often represented as a list of vertices along with a list of faces, storing vertex pointers. For the cloth simulator, this will not suffice since we will need the ability to access an edge’s neighbouring faces in a fast manner. This would be achieved by storing the topological information using a winged-edge mesh structure. To make it possible to conveniently perform additional queries like fetching a face’s adjacent faces, a half-edge mesh structure is used instead. This is a representation that stores information based on half-edges, i.e. directed edges. They are organized so that they form a linked list around the face they store a pointer to. If counter-clockwise convention is used, every half-edge stores a pointer to its left face and a pointer to the next half-edge. It also stores a pointer to its starting vertex, and a pointer to the half-edge running in parallel to it, in the opposite direction. Using this technique, it is easy to quickly perform actions like looping through all the stencils consisting of two adjacent triangles, centered around an interior edge of the mesh.
Chapter 3. Implementation of the cloth simulator 6
3.3
Backward Euler scheme
As is the case with most models for cloth simulation, the one in Baraff and Witkin’s paper is formulated as a partial differential equation which is discretized and solved as a linear system. Let x be a vector that holds the global position coefficients for all n particles, one particle at a time: x = (x0, x1, ..., xn−1)
T
∈ R3n, so that x
i ∈ R3 is the world space
position of the i’th particle. On the same format we have a net force vector f consisting of both internal and external forces acting on the cloth particles. The motion in time for a particle i is given by the equation ¨xi = fi/mi, where mi is the point mass of the
particle. The mass is calculated as the sum of the contributing mass from all triangles that contain the particle. Each of these triangles contribute to the particle’s mass with one third of the triangle’s mass, which is the cloth density multiplied by the fixed area in the planar uv-coordinate system that will be discussed further below. By defining the mass matrix M = diag (m0, m0, m0, m1, m1, m1, ..., mn−1, mn−1, mn−1), the motion in time for
the entire system is given by the ordinary differential equation M¨x = f (x, ˙x).
By defining ˙x = v, this can be written as d dt x v = v M−1f (x, v) . (3.1)
There are several possible integration schemes available that can be used on 3.1. As men-tioned above, Baraff and Witkin proposed using the backward Euler scheme, which is the approach followed here. By letting functions with subscripts n and n + 1 denote the func-tion at time t = tn and t = tn+1= tn+ h respectively, the implicit backward Euler scheme
applied to 3.1 will yield ∆x ∆v := xn+1− xn vn+1− vn = h vn+ ∆v M−1f (x n+ ∆x, vn+ ∆v) . (3.2) This looks quite similar to the explicit forward-backward Euler scheme, but has the sig-nificant difference that instead of relying on values at the current time t = tn, it also
includes values at the end of the time step, t = tn+1. 3.2 is a nonlinear equation in ∆x
and ∆v. Instead of solving this exactly, which would require iteration, Baraff and Witkin presents what some authors refer to as a semi-implicit version of 3.2, by using a Taylor series approximation to replace f (xn+ ∆x, vn+ ∆v) with fn+∂x∂f∆x +∂v∂f∆v. This results
in ∆x ∆v = h v n+ ∆v M−1(fn+∂x∂f∆x + ∂v∂f∆v) . (3.3)
Substituting ∆x with h(vn+ ∆v) in the bottom row and rearranging, we are left with
(I − hM−1∂f ∂v− h 2M−1∂f ∂x)∆v = hM −1(f n+ h ∂f ∂xvn), (3.4) where I is the 3n by 3n identity matrix. To animate the cloth, we solve for ∆v, calculate ∆x = h(vn+ ∆v) to get the state one step forward in time, and then update the mesh
with the new particle positions. In order to do this, the forces f , and the force derivatives ∂f /∂x and ∂f /∂v, have to be calculated at the beginning of each time step. (∂f /∂v is zero for the internal forces that do not depend on v, but not for the dissipative forces.)
scribe the cloth’s material behaviour is by defining a scalar-valued function over its domain, describing the potential energy, W (x). The common approach here is to separate the energy corresponding to in-plane deformation, and the energy due to out-of-plane deformation, i.e. W (x) = Wplanar(x) + Wbend(x). These measures are sometimes distinguished as
intrin-sic or membrane energy for the in-plane motions, and extrinintrin-sic or flexural energy for the out-of-plane motions. To further ease the implementation, a locally defined function, valid for only a subset of particles can be used. Then, the energy is expressed as a sum of local contributions, i.e. W (x) = P
iWi(x). The force that arises is the negative gradient of
potential with respect to position:
fi= −
∂Wi(x)
∂x .
This means that the physical modeling of the cloth is such that modifications causing it to be stretched, sheared or curved, will yield forces striving to restore the deformed state to its undeformed rest state. As shown above, the forces’ derivatives with respect to positions are also present in the system, denoted ∂x∂f. This is the force Jacobian - a square 3n × 3n matrix. Since it really consists of second order partial derivatives of the energy function, it can also be seen as the energy Hessian.
3.4.1
Potential energy and condition functions
Instead of using energy functions directly to derive the forces, Baraff and Witkin use various condition functions, C(x), formulated such that they evaluate to zero when the cloth is in its undeformed state. Their motivation for this is the difficulties of deriving sensible damping functions from energy functions. They use three different condition functions to derive the internal forces due to stretching, shearing and bending separately. To derive the forces, they define energy functions as expressions of the corresponding condition functions as
Ec(x) =
k 2C(x)
TC(x),
with k being a stiffness parameter, specific for the condition function. This gives the force contribution fi= − ∂Ec(x) ∂xi = −k∂C(x) ∂xi C(x), (3.5)
for all particles i upon which C(x) depends. For every other particle, this particular force is zero. The corresponding force derivative follows from
∂fi ∂xj = −k ∂C(x) ∂xi ∂C(x) ∂xj +∂ 2C(x) ∂xi∂xj C(x) . (3.6)
The latter two formulas correspond to equation (7) and (8) in [1]. To clarify the notation, recall that C(x) represents any of the three condition functions (corresponding to bend, stretch and shear energy). The final contributing force on a specific particle, will be the sum of the force contributions from all three separate condition functions. Important to note as well, is the fact that a particle may be affected by several different locally defined C(x), if the union of their respective domains include this particle. The three following sections will give the definition of the condition functions as given in [1].
Chapter 3. Implementation of the cloth simulator 8
3.4.2
The bend condition function
The condition function corresponding to the bend energy is defined in a local setting of two adjoining triangles, consisting of four vertices. The condition function for bending is simply defined as
Cbend(x) = θ,
with θ being the angle between the normal vectors of the two triangles (the dihedral angle). If n1 and n2 are the unit normal vectors, and e is the edge vector separating the two
triangles, θ can be computed using sin θ = (n1× n2) · e and cos θ = n1· n2. To maintain the
discrimination beteween positive and negative angles, this is implemented using the atan2 function.
3.4.3
The stretch condition function
Baraff and Witkin use planar coordinates for defining their stretch and shear conditions. This means that every particle i in the cloth should have a fixed coordinate (ui, vi) in a
planar coordinate system. Since this particular cloth simulation model is defined for planar triangle meshes only, this imposes no problem, and is implemented simply by assigning the planar coordinates to the vertices in the mesh structure. Stretch or compression is measured by looking at the derivatives of a function w(u, v) that maps from planar coordinates to world coordinates. The derivatives, wu = ∂w/∂u and wv = ∂w/∂v, have a magnitude
of 1 when the material is unstretched in the u and v directions respectively. w(u, v) is approximated to be constant over an entire triangle. To evaluate wuand wv for a triangle
with particles 0, 1 and 2, first the following quantities are evaluated: ∆x1 = x1 − x0,
∆x2= x2− x0, ∆u1= u1− u0, ∆u2 = u2− u0, ∆v1 = v1− v0 and ∆v2= v2− v0. Then
the equation wu wv = ∆x1 ∆x2 ∆u 1 ∆u2 ∆v1 ∆v2 −1 (3.7) is solved for wuand wv. This is eqn. (9) from [1]. All these quantities are readily available
for look-up in the geometric structure. Since the rightmost matrix is constant over time, because of the fixed uv coordinates, it is only computed once at the beginning of the simulation.
The condition function is finally defined as Cstretch(x) = Cu Cv = A kwu(x)k − bu kwv(x)k − bv . (3.8)
where A is the area of the triangle in uv coordinates. The parameters bu and bv can be
used to arbitrarily modify the local rest state for stretching. This is suggested to be used for modelling a desired behaviour of clothing. If for example the value of bu > 1, then the
rest length in the u direction is increased. In the current implementation, these variables are however always set to 1.
3.4.4
The shear condition function
The evaluation of the shearing in the plane of the cloth is based on the angle in world coordinates between the u and v axes of a triangle. For a triangle that is not sheared, this
Cshear(x) = Awu(x) · wv(x), (3.9)
where A is the area of the triangle in uv coordinates.
As the shear calculations for a triangle share wuand wvwith the stretch calculations, they
are computed only once in a time step, to save computation time.
3.4.5
Differentiating the condition functions
As described in 3.4.1, evaluating the force gradients involves taking the first and second order partial derivatives of the condition functions described in 3.4.2 - 3.4.4. This requires taking the second derivatives of vector valued functions with respect to vector quantities, yielding rank 3 tensors. However, the condition functions corresponding to bend and shear energy are in fact scalar valued. Only the condition function that corresponds to the stretch energy is a vector valued expression, and it can be broken down to scalar components, Cu(x) and Cv(x) during the calculations. By doing so, the calculation of ∂2C(x)/∂xi∂xj
will yield results of no higher rank than regular matrices. Although, when evaluating the derivatives analytically, the calculations are done in a component-wise fashion, in order to avoid quantities of higher order than matrices during the entire chain of expressions.
3.4.5.1 Finite differences
The initial approach is to estimate the bending related derivatives using finite differences. This is the method of approximating the derivatives of a known function, using only values of the function itself at certain discrete points. In for example [6], standard first and second order centered approximations are derived. These techniques are extended to differentiation with respect to vector quantities. The calculations are done in a component-wise fashion. For example, the formula used for first order differentiation of a scalar valued function f (x0, x1, x2, x3) with respect to component x of vertex x0is
∂f ∂x0x
= 1
h(f (x0+ hx, x1, x2, x3) − f (x0− hx, x1, x2, x3)) ,
where hx = h2( 1 0 0 )T, with h being the spatial step size. Similarly, derived as a
centered difference of centered differences, the second order partial derivative of f with respect to component x of vertex x0 and component y of vertex x1 is estimated as
∂2f ∂x0x∂x1y = 1 h2(f (x0+ hx, x1+ hy, x2, x3) − f (x0+ hx, x1− hy, x2, x3) −f (x0− hx, x1+ hy, x2, x3) + f (x0− hx, x1− hy, x2, x3)), where hy = h2( 0 1 0 )
T. The results are assembled together to yield the vector valued
forces and the matrix valued force derivatives. When the technique is employed to calculate the bending forces and bending force gradients for a local setting of two triangles meeting at a hinge edge (with four vertices involved), there are 16 different 3 × 3 matrices that need to be computed for the second order derivatives. The forces, that are the first order derivatives, result in 4 different vectors in R3.
Chapter 3. Implementation of the cloth simulator 10 3.4.5.2 Analytic expressions
Dean Macri at Intel wrote a paper [7] in which he derived expressions for the analytic derivatives of all condition functions from Baraff and Witkin’s paper. Some minor mistakes in his formulas were pointed out in an unpublished article by David Pritchard [8]. With these corrected, Macri’s derived equations could be used to implement the calculations. However, they necessitate the use of rank 3 tensors for some of the derivatives. To ease the implementation, the calculations can instead be done for one component at a time. This is the approach used in [8], which contains a detailed derivation of the analytic derivatives of the condition functions, that is being consulted in this project’s implementation. However, since the aim of this project is to fully replace the bending calculations, the finite difference approach is considered good enough. Only the stretch and shear related calculations are implemented as analytical expressions. The derivatives used in the implementation will be presented below.
Following from equation (3.7) are the quantities wu= (x1− x0)∆v2− (x2− x0)∆v1 ∆u1∆v2− ∆u2∆v1 , wv= −(x1− x0)∆u2+ (x2− x0)∆u1 ∆u1∆v2− ∆u2∆v1 .
In order to take the derivatives of the stretch condition function in equation (3.8), the following derivatives are calculated for wu
∂wu ∂xi = ∂wux ∂xix I, ∂wux ∂x0x = ∆v1− ∆v2 ∆u1∆v2− ∆u2∆v1 , ∂wux ∂x1x = ∆v2 ∆u1∆v2− ∆u2∆v1 , ∂wux ∂x2x = −∆v1 ∆u1∆v2− ∆u2∆v1 . The corresponding expressions for wv are
∂wv ∂xi =∂wv x ∂xix I, ∂wv x ∂x0x = ∆u2− ∆u1 ∆u1∆v2− ∆u2∆v1 , ∂wv x ∂x1x = −∆u2 ∆u1∆v2− ∆u2∆v1 , ∂wv x ∂x2x = ∆u1 ∆u1∆v2− ∆u2∆v1 .
∂xi = A ∂xix wu, ∂Cv ∂xi = A∂wv x ∂xixwv, ∂2C u ∂xi∂xj = A kwuk ∂wux ∂xix ∂wux ∂xj x(I − wuw T u), ∂2C v ∂xi∂xj = A kwvk ∂wv x ∂xix ∂wv x ∂xj x(I − wvw T v).
When computing equation (3.6) from above, instead of directly performing the multiplica-tion ∂2Cstretch(x)
∂xi∂xj Cstretch(x), which would mean multiplying a 3n × 3n × 3n quantity with a 3n vector, it is done one component at a time:
∂2C stretch(x) ∂xi∂xj Cstretch(x) = ∂2C u(x) ∂xi∂xj Cu(x) + ∂2C v(x) ∂xi∂xj Cv(x)
The derivatives of the shear condition from equation (3.9) function use the above derivatives of w as well, and are calculated as
∂Cshear ∂xis = A(∂wu ∂xis · wv+ wu· ∂wv ∂xis ) = A(∂wux ∂xix wv s+ wus∂wv x ∂xix ), ∂2Cshear ∂xis∂xj t = A( ∂ 2w u ∂xis∂xj t · wv+ ∂wu ∂xis ·∂wv ∂xj t+ ∂wu ∂xj t · ∂wv ∂xis + wu· ∂2wv ∂xis∂xj t ) = ( A(∂wu x ∂xi x ∂wv x ∂xj x + ∂wu x ∂xj x ∂wv x ∂xi x) : s = t, 0 : s 6= t.
where s and t are vertex components {x, y, z}. Note that the second order derivatives of w are zero.
3.5
Dissipative forces
The damping forces that are introduced by Baraff and Witkin are, like the condition func-tions’ resulting forces, acting only in the direction of the gradient of the condition functions. They let the damping forces depend on the system’s velocity in that direction. The reader is referred to the original article for the full details. After leaving out some terms in their equations that break the symmetry needed for insertion into the system of equations, and simplifying, the following are inserted into the global versions in equation (3.4):
fdi= −kd ∂C(x) ∂xi ˙ C(x), ∂fdi ∂xj = −kd ∂2C(x) ∂xi∂xj ˙ C(x), ∂fdi ∂vj = −kd ∂C(x) ∂xi ∂C(x) ∂xj T ,
Chapter 3. Implementation of the cloth simulator 12 where the time derivative of C(x) is evaluated as
˙ C(x) =X i ∂C(x) ∂xi T ˙xi.
This requires only adding and multiplication of quantities already calculated.
3.6
Solving the system of linear equations
Recall the formulation of the system from equation (3.4): (I − hM−1∂f ∂v− h 2 M−1∂f ∂x)∆v = hM −1(f n+ h ∂f ∂xvn).
When solving for ∆v it is done in an iterative way, generating a series of approximations of a solution. The matrix on the left hand side (to the left of ∆v) is involved in one or several matrix-vector multiplications in each of the iterations. The number of iterations required to reach a solution within a given desired range of accuracy, will depend on properties of the matrix.
There are several methods available for this iterative solver approach. If trying to solve this using a regular Conjugate Gradient (CG) method, it has to be considered that these algorithms need the matrix to be symmetric, and preferably positive definite as well. This can be achieved by multiplying both sides of the equation with the mass matrix M, which results in the linear system
(M − h∂f ∂v− h 2∂f ∂x)∆v = h(fn+ h ∂f ∂xvn), (3.10) having the same solution as equation (3.4). However, as will be described in the next section, the article of Baraff and Witkin is altering the equation when introducing techniques for imposing constraints on specific particles. This results in a system that cannot easily be made positive definite. Therefore, they ultimately choose not to use a regular CG method on their modified version of (3.4). Instead they propose a modified version of the CG method, which can operate on (3.10), and in which they have treatment of constraints built in. The constraint-specific data is provided to the solver as an additional matrix and a vector.
3.6.1
Constraints
To be able to incorporate the response of collisions into the simulation, Baraff and Witkin propose different techniques for the self-collisions between different parts of the same cloth, and the collisions between the cloth and solid objects in the same environment. For self-collisions they use strong springs and position alteration, and for cloth/solid self-collisions they use what they call mass modification. The self-collision response is completely omitted in this thesis’ final implementation, as described in section Limitations. For imposing contact constraints on individual particles, as a result of collisions with solid objects, the techniques of mass modification along with position alteration is implemented.
Mass modification uses the idea that the mass of a particle is altered so that the resulting acceleration can be controlled as desired. Of course the strict interpretation of mass is a scalar, but we have in the model a mass matrix with 3 components for every particle. This means that by altering different components of the mass matrix, we would be able
¨ xi= 0 0 0 0 1/mi 0 0 0 1/mi fi. (3.11)
The matrix can be seen as the block from M−1that corresponds to particle i. By having set the first element in the block to zero, we force the acceleration in the x direction to be zero. We are not limited to just locking a specific axis, but can rather decide the acceleration in an arbitrary direction. If the inverse mass matrix in the expression (3.11) above would have remained unaltered, it would be equal to m1
iI, where I denotes the 3x3 identity matrix. If we instead use an inverse mass matrix of m1
i I − pp
T, with p = 1 0 0 T
, we would of course get the same result as equation (3.11). This modified inverse mass matrix is an orthogonal projection matrix, and by setting these individually for the particles, we have a way of controlling the motions as intended.
To generalize this, Baraff and Witkin define the system’s full modified version of the inverse mass matrix M−1 as W. By letting S being a matrix with its diagonal blocks Si being the
desired orthogonal projection for particle i, i.e. S = diag(S0, ..., Sn−1), the corresponding
diagonal blocks of W are given by Wi=m1
iSi.
The different possibilities for projections are broken down into the following cases:
Si= I : ndof (i) = 3, (I − pipTi) : ndof (i) = 2, (I − pipTi − qiqTi) : ndof (i) = 1, 0 : ndof (i) = 0,
where ndof (i) means particle i’s number of degrees of freedom. The vectors p and q are the constrained directions. They are mutually orthogonal unit vectors. By defining a vector z that contains the desired solution in the constrained directions, the velocity change is not limited to be zero for the constrained particles. To clarify, this would conceptually mean that we switch to using the system
(I − hW∂f ∂v− h 2W∂f ∂x)∆v = hW(fn+ h ∂f ∂xvn) + z. (3.12) Note that the vector z needs components equal to zero in the unconstrained directions. However, as mentioned in the previous section, the matrix here cannot easily be made positive definite by multiplication with the mass matrix. Therefore, Baraff and Witkin ultimately choose not to use a regular CG method on a modified version of (3.12). Instead they propose a modified version of the CG method, which can operate on (3.10), and in which they have the treatment of constraints built in. The constraint-specific data is provided to the solver as the matrix S and the vector z.
By formulating the problem as
SA∆v = b, (I − S)∆v = z,
Chapter 3. Implementation of the cloth simulator 14 we will have a cofiguration in which the equation of motion has an effect in the subspace projected by S. In the orthogonal subspace projected by (I − S), the vector z will give the solution.
By designing the modified CG method to never allow an update in the constrained direction, it will reach a solution in which the constrained directions will be equal to the components in the vector z. The article calls the proposed algorithm the Modified Pre-conditioned Con-jugate Gradient (MPCG) method. In this thesis implementation, another version described in the next section is used.
3.6.2
Ascher and Boxerman’s modified pre-conditioned conjugate
gradient method
In 2003, Ascher and Boxerman proposed an improved version of the MPCG method with reportedly better performance [5]. They also provide a proof of convergence. This is the algorithm used in the implementation and it is reproduced in detail here.
With A = M − h∂f ∂v− h
2 ∂f
∂x and b = h(f0+ h ∂f
∂xv0), from equation (3.10), the linear system
can be written as A∆v = b. We need to solve for ∆v, and keep the soution to obey the constraints imposed on the particles. This information is available to the solver via the filtering matrix S, and the vector z, which contains the desired solution for the constrained directions.
The algorithm is as follows:
Provided the known quantities A, b, S and z are assigned their proper values, we start by setting the initial iterate to ∆v = Sy + (I − S)z. The vector y is a guess for the case of no constraints. This can be an all zero vector, or as chosen in this implementation, the result from the last time step. Next we initialize the following variables:
ˆ b = S(b − A(I − S)z), bδ = bˆTP−1ˆb, r = S(b − A∆v)), p = SP−1r, δ = rTp. Then, while δ > tolerance2b
δ, we do: s = SAp, α = δ/(pTs), ∆v = ∆v + αp, r = r − αs, h = P−1r, δold = δ, δ = rTh, p = S(h+(δ/δold)p).
The number of iterations required to reach a solution within the desired interval varies depending on properties of the matrix A. The choice of pre-conditioner matrix should affect the performance as well. There are several different approaches in the literature regarding
Functionality for setting up the filtering matrix and the constraint vector is implemented. Figure 3.2 below is showing a screen capture from a simulation where the particles along one of the border edges of a piece of cloth are completely locked. Even though external forces (i.e. gravity and the simple wind model described above/below) are present, the solution for the locked particles remains zero.
Figure 3.1: Snapshot from a simulation in which a piece of cloth is completely constrained along one border edge. Gravity and wind are present.
3.6.3
Exploiting the matrix sparsity
The discretization of the physical system will lead to a matrix representation with a low proportion of nonzero elements. This is something that can be exploited to save storage space and be able perform more efficient calculations. For very large matrices, taking advantage of the sparsity may be necessary to be able to store the matrices in memory. There are different techniques for storing a sparse matrix. This implementation uses UBLAS from the Boost libraries. Using a thoroughly tested existing library for the sparse matrix and vector operations eases debugging of the code. Due to the slight overhead costs of UBLAS, there are probably some performance advantages to be gained for the solver if it is rewritten with more application specific expressions. There are several other free linear algebra packages available to choose from as well.
Chapter 3. Implementation of the cloth simulator 16
3.7
Collision handling
The handling of collisions for the cloth can be divided into two quite different parts; self-collisions and self-collisions with solid objects. A cloth where some parts are colliding with other parts of the same cloth is generally much more difficult to handle. Contrary to collisions with solid objects, we cannot simply check if a particle is inside or outside the other object, since the cloth is an infinitesimally thin surface. The response code is harder to construct as well. If we have a large chunk of wrinkled cloth involved in a large amount of collisions, then an impulse aimed to resolve one collision may very well lead to new collisions. Collision response for contact with solid objects is implemented in the simulator using the technique for imposing constraints as described above, together with position alteration to resolve any collided particles. This is how Baraff and Witkin propose to move a particle as desired. If we were to simply change the position of a particle after a step has been taken, it would result in strange artefacts of jumping particles. The neighbouring particles would receive no information in advance, so when moving the particle it would lead to too large stretch energies. If the adjustment vector for the particle positions, called y, is instead incorporated into the equation for the change in positions as
∆x = h(vn+ ∆v) + y,
then equation (3.10) would instead look like (M − h∂f ∂v− h 2∂f ∂x)∆v = h(fn+ h ∂f ∂xvn+ ∂f ∂xy). (3.13) By using this technique, a particle can safely be moved to resolve a collision in one time step. For rendering, it is best to move all vertices to a collision free state in a post-processing step that does not affect the physical position vector. Below is a sequence of images showing the collision with a solid sphere.
3.8
The simulation loop
Described above are the quantities being calculated each time step of the simulation. The sparse matrix and dense vector from equation (3.13) are assembled by adding the force and force Jacobian contributions from one type of condition function at a time, by looping over the local domains. External forces and gradients are added as well. The matrix S, which contains the constraints that are to be imposed on certain particles, is assembled after collision detection. So is the vector that contains the values for position alteration of collided particles. When the system is set up, the solver iterates until a solution for ∆v within the desired accuracy is reached. Then the new positions are updated. Any new collisions occurring are handled by non-physically moving such particles, after which the frame can be rendered.
3.9
Conclusion
Using finite differences to approximate the derivatives of the bend condition function in order to be able to evaluate the forces and the force gradients, proved to be accurate enough to get a working simulation up and running. However, it is not very effcient due to
Figure 3.2: Snapshots from a simulation. A piece of cloth, constrained along one edge, falling over a solid immobile sphere. The mesh consists of 200 triangles and 121 vertices. The duration of the sequence is 55 seconds, with a frame rate of 8-10 fps on a 2.2 GHz Intel Core 2 Duo CPU.
the large number of operations required at each time step. Finding the optimal step size is not so obvious either. One benefit of this approach is the generality that it offers. When experimenting with different models, it may be convenient to not have to derive an analytic expression for the derivatives of every function. This applies, of course, especially when the function is non-linear in the independent variables which we are to differentiating with respect to, and thus complex to differentiate.
When simulating a cloth consisting of 121 particles as in figure 3.2, the computation time is about 0.11 seconds per frame, using a 2.2 GHz Intel Core 2 Duo CPU. Introducing an-alytical derivatives for all forces would probably yield better performance, both regarding computation time and stability. During simulations with the implemented numerical ap-proximations, some unexpected behaviour in the simulation has been noticed on occasion. However, the implementation of the simulator described above yielded a fairly robust envi-ronment for trying out different alterations of the models for both the bend resisting forces, as well as the in-plane stretch and shear resisting forces. The simulation remains stable, even with relatively large time steps. Although, to guarantee stability, making the time step adaptive would be a good precaution. This can be done by increasing the step size temporarily whenever the planar strain in the cloth exceeds a certain amount. After a number of subsequent successful steps, the step size is gradually restored (see e.g. [4]).
Chapter 4
Modification of the bending
model
This chapter will describe the substitution of the current bending model from [1], with two more recently proposed models. The first model is fully derived in the article Discrete quadratic curvature energies [9] and is suited for surfaces with a flat rest state, e.g. cloth and thin plates. The second model, presented in Cubic Shells [10], is capable of simulating curved surfaces as well. This opens up the possibility of simulating thin shell objects, such as bowls, hats, leaves etc.
4.1
Quadratic Curvature Energies
In [9], Bergou et al presented a bending model derived by defining bending energy as a curvature measure. In that sense it is not very different from earlier models. What is new, is the recognition of the energy being a quadratic polynomial in positions, for the special case of isometric surface deformations (where the surface is allowed to bend but not stretch or shear). They relax the demand of isometry, and state that their model is valid for inextensible plates where the membrane stiffness is greater than the bending stiffness by more than four orders of magnitude. To derive their discrete bending model, they start out in the continuous setting, defining the bending energy of a deformable surface as an integral of the squared mean curvature over the surface area:
Eb(S) = 1 2 Z S H2dA. (4.1)
Important properties of the energy formulation are the invariance under rigid motions (since the internal energy should not depend on the embedding in space), and uniform scaling of the surface.
If the embedding map of the surface S in R3
is written as x : S → R3, the mean curvature
normal of S can be written as the Laplace-Beltrami operator ∆ applied to the embedding: H = ∆x.
S
The article refers to this equation, in the context of isometric deformation, as the isometric bending model (IBM). They proceed by defining a discrete version, which shares the prop-erty of being quadratic in positions. This implies a linear gradient, and a constant Hessian, which obviously means a significant decrease in the number of computations required each time step of a simulation.
With the same notation for the vertex position vector x as in previous sections of this thesis, the article’s discrete version of Eb is written as
Eb(x) = 1 2x TQ x = 1 2 X Qijhxi, xji.
Q is a quadratic form that has to maintain the properties discussed above. The matrix is defined as
Q = LTM−1L,
where L is the discrete Laplacian and M is a mass matrix. By noting that M−1Lx is a
discrete pointwise version of the mean curvature vector ∆x, and that the mass matrix can be interpreted as a discrete ”area measure”, we get the expression for Q written above by looking at the analogous discrete version of 4.2
Eb(x) = 1 2(M −1Lx)TM(M−1Lx) = 1 2x T(LTM−1L) x.
The reader is referred to the original article for a thorough discussion of how they choose a suitable discretization of the Laplacian and the mass matrix. For the application of cloth and plate simulation, they ultimately propose a version of the discrete Laplacian centred around the edges. (When the variables in calculations on a discretized surface are placed at the edge midpoints rather than at the vertices, it is referred to as using non-conforming elements - a technique employed e.g. in Finite Element Method calculations.) They then use a mapping from vertices to edge midpoints, in order to combine the edge-based Laplacian with vertex-edge-based degrees of freedom, since the resulting forces and force gradients are required at each vertex in order to use it in the regular simulation context with a continuous mesh grid.
4.1.1
Implementation
Even though the original article provides all the necessary formulas to compute the Hessian matrix Q, in an accompanying paper [11] they provide a slightly different but equivalent approach. The mapping, Laplacian and mass matrix are simplified to a nice expression for the 4 × 4 sub matrix Q(e0), local for an internal edge. The computations are being done
for a two-triangle subspace, like what was done with the bend calculations in chapter 3. ¯ Qi= 3 ¯ Ai ¯ KTiK¯i. ¯
K is the row vector (¯c03+ ¯c04, ¯c01+ ¯c02, −¯c01− ¯c03, −¯c02− ¯c04), with ¯cjk = cot ∠¯ej, ¯ek.
The local matrices are assembled together to yield the global bend energy Hessian. The resulting force contributions are computed as fbend= Qx. The matrix Q is in Rn× Rn, so
Chapter 4. Modification of the bending model 20 every third element), or as done in this implementation, constructing a R3n× R3n version
of Q with each element repeated three times. This way it is easy to insert the bending contribution into the global force Jacobian (∂f∂x), by adding the two matrices together. When using this bending model instead of the previous one, the bend energy Hessian (the 3n × 3n version) is pre-computed once, before the simulation starts. The global matrix ∂x∂f is updated with these values. Then, for each time step, the bending forces are computed by the global matrix-vector multiplication. The time integration technique, as well as the model for in-plane motions are kept the same as before.
4.1.1.1 Results
By replacing the bending calculations of Baraff & Witkin with this new model, the most noticeable effect is the drastic decrease in computation time at each simulation time step. Since the force gradient (energy Hessian) computations were the most time consuming of the previous version of the simulator, it is a great relief to have a constant Hessian. The force computations are faster as well. The simulation is stable, and compared to the previous model, it is now possible to have stiffer bending forces than before. Air drag forces prove to be enough when it comes to damping, even though the article suggests using dissipative forces fd= Qv. Figure 4.1 below is showing a sequence of images where a square patch of
thin material is colliding with a sphere. Relatively large stiffness coefficients make it rather reluctant to bend. Noticeable from the figure is the fact that the discretization of the surface has quite an impact on how close the behaviour in the continuous setting is actually captured. For best results, the edge directions in the plane would be evenly distributed.
Figure 4.1: Snapshots from a simulation, spanning over 10 seconds. A thin patch of a semi-stiff material is falling down over a solid sphere. The surface will seek its rest state after impact. The mesh consists of 200 triangles and 121 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 26-28 fps.
curvature is not equal to zero), a recent article called Cubic Shells arrived at a bend energy formulation not very different from the IBM described in the previous section [10]. They extend equation (4.1) above to instead use the squared change in mean curvature from undeformed to deformed state:
Eb(x) = 1 2 Z S (H − ¯H)2dA.
A bar over a quantity will here and henceforth indicate that it is to be evaluated in the undeformed setting. Whereas the article from the previous section rearranged the energy formulation to a quadratic polynomial in x under isometric deformations, the corresponding version in [10] is rewritten as Eb(x) = 1 2 Z S hH, Hi − 2H, ¯H ˆn + ¯H2 dA.
They show that this is a cubic polynomial in x, which will lead to corresponding forces that are quadratic. They then go on to derive a discrete version, with this property in common. The discretization of the energy is carried out by looking at the hinge based energy
Eb(x) = 1 2 X i 3k¯eik2 ¯ Ai 2 sinθi− ¯θi 2 2 ,
with the sum being taken over all interior edges i. The contribution from each interior edge can be rearranged and written as
Eb(x)i=
3k¯eik2
¯ Ai
1 − cos θi− ¯θi , (4.3)
which in fact is a cubic polynomial in x, even though as pointed out in the article, the above energy coincide with the previous models from the literature (using 2 sin(θ/2), 2 tan(θ/2) or θ2 as in chapter 3), when θ approaches 0 as it does when performing a proper mesh refinement.
The article introduces the vectors t(0)0 and t(1)0 , perpendicular to, and of the same length as e0 and e1. They can be computed as
t(0)0 = − cot α03e1− cot α01e3,
t(1)0 = − cot α04e2− cot α02e4.
(4.4)
Defining β = ke1
0k(cot α01+cot α03)(cot α02+cot α04), equation (4.3) is equal to the following expression, after some rearranging and employing of the cosine identity
3 ¯ A0 k¯e0k2+ ht (0) 0 , t (1) 0 i cos ¯θ +3 ¯¯β A0 (e0× e1) · e2 sin ¯θ.
The first term corresponds to the thin plate energy, and is equivalent to the energy from the previous section. The second term is called the cubic term, and vanishes for flat rest states.
As in the article with the previous model, a simplified expression for calculating the forces and force gradients resulting from this energy is given. This model’s force Jacobian can be
Chapter 4. Modification of the bending model 22 pre-computed once as well.
4.2.1
Implementation
Quite similar to the previous section, for every interior edge i, a local 4 × 4 matrix ¯Qi is
calculated before simulation starts as ¯ Qi= 3¯λicos ¯θi ¯ Ai ¯ KTiK¯i.
The row vector ¯K is set up exactly as before: ¯K = (¯c03+¯c04, ¯c01+¯c02, −¯c01−¯c03, −¯c02−¯c04),
with ¯cjk = cot ∠¯ej, ¯ek. ¯λi is a stiffness value that can vary among the hinge edges. Since
no anisotropic bending is used in this implementation, they are all set to the same value. The global n × n matrix ¯Q is assembled in the same way as before, and contributes to the global force Jacobian every time step. The forces arising from this model are separated into thin plate forces and non-flat forces. Since the only difference of this energy Hessian and that of the previous model, is the scaling factor cos ¯θ that is maximised for flat rest states, the thin plate forces are calculated as previously. The non-flat forces are calculated as
f0= −f1− f2− f3,
f1= ¯k(e1× e2),
f2= ¯k(e2× e0),
f3= ¯k(e0× e1),
with the coefficient ¯
k = 3λi(¯c01− ¯c03)(¯c04− ¯c02)((e0× e1) · e2)/( ¯Aik¯eik3). (4.5)
For a better understanding of how these forces act together, see figure 4.2 below for an illustration of their interplay in a local setting of one hinge edge. The blue lines represent the forces striving to achieve a flat state, the red lines represent the counteracting forces that, on their own, would yield an angle of π/2. Together, they give rise to the resultant forces, represented by the bold green lines, that will act to achieve the non-flat rest state. Note that this figure is merely produced for helping the conceptual understanding, and the scaling and direction may not exactly match the result as given by the unaltered quotients of the force coefficients from the original article. Note also that the model is not necessarily valid for local deformations as large as those displayed in the figure.
4.2.1.1 Changing in-plane energy model
In order to simulate naturally curved surfaces, the model for in-plane motions needs to be updated as well. If every vertex is assigned additional sets of uv-coordinates, so that different ones are used depending on which triangle is evaluated, the first model could still be used. In this implementation however, a system of springs is used with each edge in the mesh acting as a linear spring between its two endpoint vertices. The energy function for the spring model is borrowed from Choi and Ko’s article Stable but responsive cloth [12]. Defining xij = xi− xj and ¯xij = ¯xi− ¯xj , the formula is written
E = (1 2ks(kxijk − k¯xijk) 2 : kxijk ≥ k¯xijk, 0 : kxijk < k¯xijk.
Figure 4.2: Visualization of the forces resulting from the cubic bending model, acting on the vertices in a two-triangle hinge, having a rest angle of π/9. The angles have been deformed to π/4 and 3π/4 radians respectively, for the left and right image. The non-flat rest state is shown by the dashed green lines, and the flat state is shown by the dashed blue lines. The solid blue lines represent the forces acting to achieve a flat state, whereas the solid red lines depict the counter-acting non-flat forces. The resultant forces, shown by solid green lines, will act to achieve the non-flat rest state.
The force acting on particle i and the force Jacobian are computed as fi= ∂E ∂xi = ( ks(kxijk − k¯xijk)( xij k¯xijk) : kxijk ≥ k¯xijk, 0 : kxijk < k¯xijk. ∂fi ∂xj = ks xijxTij xT ijxij + ks(1 − k¯xijk kxijk)(I − xijxTij xT ijxij) : kxijk ≥ k¯xijk, 0 : kxijk < k¯xijk. 4.2.1.2 Adding plasticity
In the real world, there exist no materials that are completely elastic. Of course, for many animation applications, an elastic model can produce animations that are good enough, even if it is not very true to nature. But there are numerous cases where the plastic behaviour of a material is crucial for the realism of an animation. Consider for example a ball of clay being dropped onto a hard floor, a piece of paper being crumpled or folded, or the body of a car being involved in an accident. These are situations which cause the material to be permanently deformed.
A simple model for plastic deformations is implemented, building from the techniques de-scribed in [13] and [14], where plasticity is dede-scribed in terms of strain only. Strain is a geometric measure of the deformation away from the material’s rest state, and not to be confused with stress, which describes the amount of load exerted on the material. The total strain of an object can be seen as the sum of the plastic strain p and the elastic strain e. When an object is subject to deforming forces, for some interval of resulting strain it
will generally deform in the so called elastic regime - when the forces are removed, it will be restored to the original rest state. For a plastic material, in this rest state the strain is equal to the plastic strain. If deforming forces are increased so that the strain exceeds a certain threshold, any further increase in total strain will add to the plastic strain. This means that when the forces that have resulted in a plastic deformation are removed, the
Chapter 4. Modification of the bending model 24 object will seek a new rest state different from the original one. We generally do not want to update the plastic strain with an amount equal to the strain above the threshold, since a realistic behaviour would allow the material to be ”hardened”. Instead, as described in [14], the plastic threshold is set to grow linearly with the accumulated plastic strain: plastic regime starts when the elastic strain e = − p reaches 0+ γp, where 0 is the initial
threshold value and γ is a constant. The change in plastic strain is expressed as
∆p = ( − 0)/(1 + γ). (4.6)
In the implementation, the following steps are taken each time the bending forces are evaluated for a hinge edge i: The elastic strain is evaluated as e = k¯Aeik
i (θi− ¯θi), using the same notation as introduced in previous sections. The plastic strain is expressed as p = k¯eik
Ai (¯θi). With 0 and γ being user defined parameters, the threshold is evaluated as 0+γp. If the elastic strain exceeds the threshold, the plastic strain is updated using formula
(4.6). The rest angle ¯θi for the hinge is updated as well, by setting it to be equal to pk¯Aei
ik. Since the rest state has changed, a recalculation of the energy Hessian is necessary. This is done by subtracting the old values from the corresponding components of the Hessian matrix, and adding those resulting from the new rest angle. The same thing is done for the non-flat force coefficients given by formula (4.5).
4.2.1.3 Results
Incorporated into the simulation framework, the cubic bending model provides an interest-ing way to animate curved surfaces. This change drastically increases the possible areas of application for this type of simulator. Because of the straightforward way of assembling the bending force Jacobian, the cubic model is suitable for use in any simulator that uses an implicit time integration. As the bending force Jacobian is pre-computed once, it is significantly faster than e.g. the model from [1].
Figure 4.3 is showing a sequence of renderings from a simulation without plasticity. The simulated thin shell mesh is drawn with flat shading so that the polygons are distinguishable. For a mesh of a size like this, the simulation runs at interactive frame rates. However, for simulations of larger models, real-time performance is not to be expected.
To capture a true behaviour of plastic materials, a unified handling of plasticity for both planar and bending deformations would be required. However, as illustrated by figure 4.4, the more naive way of using a fully plastic behaviour for in-plane deformation can be used to achieve the desired effect, in the case of a relaxed demand on isometry.
Figure 4.3: Snapshots from a simulation, spanning over 1.5 seconds. A thin shell object, shaped like a bowl, is falling to a solid half-sphere, fixed to the ground. The surface will seek its curved rest state after impact. The mesh consists of 130 triangles and 77 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 62-64 fps.
Chapter 4. Modification of the bending model 26
Figure 4.4: Snapshots from a simulation with plasticity, spanning over 0.2 seconds. A thin shell object, shaped like a bowl, is falling to a solid half-sphere, fixed to the ground. Notice how the dent is not fully restored after impact. The mesh consists of 130 triangles and 77 vertices. Running on a 2.2 GHz Intel Core 2 Duo CPU, the frame rate is 62-64 fps.
’98: Proceedings of the 25th annual conference on Computer graphics and interactive techniques, pages 43–54, New York, NY, USA, 1998. ACM. ISBN 0-89791-999-8. doi: http://doi.acm.org/10.1145/280814.280821.
[2] Robert Bridson, Ronald Fedkiw, and John Anderson. Robust treatment of collisions, contact and friction for cloth animation. In SIGGRAPH ’02: Proceedings of the 29th annual conference on Computer graphics and interactive techniques, pages 594–603, New York, NY, USA, 2002. ACM. ISBN 1-58113-521-1. doi: http://doi.acm.org/10. 1145/566570.566623.
[3] Matthias Teschner, Bruno Heidelberger, Matthias M¨uller, Danat Pomeranets, and Markus Gross. Optimized spatial hashing for collision detection of deformable objects. pages 47–54, 2003.
[4] R. Bridson, S. Marino, and R. Fedkiw. Simulation of clothing with folds and wrin-kles. In SCA ’03: Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 28–36, Aire-la-Ville, Switzerland, Switzerland, 2003. Eurographics Association. ISBN 1-58113-659-5.
[5] Uri M. Ascher and Eddy Boxerman. On the modified conjugate gradient method in cloth simulation. The Visual Computer, 19(7):526–531, 2003. URL http://dx.doi. org/10.1007/s00371-003-0220-4.
[6] Randall J. LeVeque. Finite difference methods for differential equations. Draft version for use in the courses AMath 585-6, University of Washington, 1998.
[7] Dean P. Macri. Real-time cloth. Intel Corporation.
[8] David Pritchard. Implementing Baraff & Witkin’s cloth simulation. May 2002. [9] Miklos Bergou, Max Wardetzky, David Harmon, Denis Zorin, and Eitan Grinspun.
Discrete quadratic curvature energies. In SIGGRAPH ’06: ACM SIGGRAPH 2006 Courses, pages 20–29, New York, NY, USA, 2006. ACM. ISBN 1-59593-364-6. doi: http://doi.acm.org/10.1145/1185657.1185663.
[10] Akash Garg, Eitan Grinspun, Max Wardetzky, and Denis Zorin. Cubic shells. In SCA ’07: Proceedings of the 2007 ACM SIGGRAPH/Eurographics symposium on Computer animation, pages 91–98, Aire-la-Ville, Switzerland, Switzerland, 2007. Eurographics Association. ISBN 978-1-59593-624-4.
[11] Miklos Bergou, Max Wardetzky, David Harmon, Denis Zorin, and Eitan Grinspun. A quadratic bending model for inextensible surfaces. In SGP ’06: Proceedings of the fourth Eurographics symposium on Geometry processing, pages 227–230, Aire-la-Ville, Switzerland, Switzerland, 2006. Eurographics Association. ISBN 30905673-36-3.
Bibliography 28 [12] Kwang-Jin Choi and Hyeong-Seok Ko. Stable but responsive cloth. ACM Trans. Graph., 21(3):604–611, 2002. ISSN 0730-0301. doi: http://doi.acm.org/10.1145/ 566654.566624.
[13] James F. O’Brien, Adam W. Bargteil, and Jessica K. Hodgins. Graphical modeling and animation of ductile fracture. ACM Trans. Graph., 21(3):291–294, 2002. ISSN 0730-0301. doi: http://doi.acm.org/10.1145/566654.566579.
[14] Yotam Gingold, Adrian Secord, Jefferson Y. Han, Eitan Grinspun, and Denis Zorin. A Discrete Model for Inelastic Deformation of Thin Shells. Technical report, 2004.
