Notes 1D Elastic Continuum From last class: elastic rod Required - - PowerPoint PPT Presentation

SLIDE 1

1 cs533d-term1-2005

Notes

Required reading:

• Baraff & Witkin, “Large steps in cloth

animation”, SIGGRAPH98

• Grinspun et al., “Discrete shells”, SCA03
• Bridson et al., “Simulation of clothing with

folds and wrinkles”, SCA03

2 cs533d-term1-2005

1D Elastic Continuum

From last class: elastic rod

• linear density “rho” (not necessarily constant)
• Youngs modulus E (not necessarily constant)
• Paratemerized by p

If homogenous, simplifies to:

• x(p) = 1
• p E(p)
• p x(p) 1
• 2x

t 2 = E

• 2x

p2

3 cs533d-term1-2005

Sound waves

Try solution x(p,t)=x0(p-ct) And x(p,t)=x0(p+ct) So speed of “sound” in rod is Courant-Friedrichs-Levy (CFL) condition:

• Numerical methods only will work if information

transmitted numerically at least as fast as in reality (here: the speed of sound)

• Usually the same as stability limit for good explicit

methods [what are the eigenvalues here]

• Implicit methods transmit information infinitely fast

E

• 4

cs533d-term1-2005

Why?

Are sound waves important?

• Visually? Usually not

However, since speed of sound is a material

property, it can help us get to higher dimensions

Speed of sound in terms of one spring (using

linear density m/L) is

So in higher dimensions, just pick k so that c is

constant

• m is mass around spring [triangles, tets]
• Optional reading: van Gelder

c = kL m

SLIDE 2

5 cs533d-term1-2005

Potential energy

Another way to look at the elastic spring forces:

through potential energy

Recall for a system at position(s) x, potential

energy W(x) gives the force

For example, this gives conservation of total

energy K+E: F = W x

d dt

1 2 vT Mv + W

( ) = vT M dv

dt + W x dx dt = vT Ma FTv = vTF vTF = 0

6 cs533d-term1-2005

Spring potential

For a single spring,

• Note were squaring the percent deformation (so this

always increases as we move away from undeformed), and scaling by the strength of spring and by the length (amount of material it represents)

To get the force on i, differentiate w.r.t. xi:

Wij = 1

2 kij

xi x j Lij 1

• 2

Lij fij = Wij xi = kij xi x j Lij 1

• xi x j

xi x j

7 cs533d-term1-2005

1D Continuum potential

Add up the potential energies for each spring to

approximate the total potential energy for the elastic rod:

Take the limit as p goes to zero: Now: how do we get forces out of this? The

negative gradient of W w.r.t. x(p)?

W

1 2 Ei+ 12

xi+1 xi pi+1 pi 1

• 2

(pi+1 pi)

i

• W =

1 2 E(p) x

p 1

• 2

dp

1

• 8

cs533d-term1-2005

Directional derivatives (regular calculus)

Pick a direction, or test vector, q Direction derivative along q is: Or alternatively And the gradient W/x is the vector s.t.

DqW (x) = lim W (x + q) W (x)

• DqW (x) =
• g (0)

where g() = W (x + q)

DqW (x) = W x iq q

SLIDE 3

9 cs533d-term1-2005

The variational derivative

We want to take the “gradient” of a continuum

potential energy to get the force density:

But how do you differentiate w.r.t. a function

x(p)?

Lets first look at a directional derivative: look at

the energy at x+q

• q is a direction, or test function

f = W[x] x where W[x] =

1 2 E x

p 1

• 2

dp

1

• g() = W[x + q] =

1 2 E x

p + q p 1

• 2

dp

1

• 10

cs533d-term1-2005

The variational derivative 2

Now differentiate w.r.t. the scalar: And evaluate at 0 to get the directional derivative

• g () = d

d

1 2 E x

p + q p 1

• 2

dp

1

• =

1 2 E d

d x p 1

• 2

+ 2 x p 1

• q

p + 2 q p

2

• dp

1

• =

E x p 1

• q

p + q p

2

• dp

1

• DqW[x] =
• g (0) =

E x p 1

• q

p dp

1

• 11

cs533d-term1-2005

The variational derivative 3

We want to make it look like an inner-product of

the “gradient” with q()

• Use integration by parts:

Ignoring boundary conditions for now, we see

that the variational derivative of W at some interior point p is just: DqW[x] =

• p E x

p 1

• qdp

1

• + E x

p 1

• q
• 1

p E x p 1

• 12

cs533d-term1-2005

Force density

The elastic force density is the negative of the

variational derivative at a point:

The acceleration of that point is force density

divided by mass density:

Which is exactly what we got before!

f = p E x p 1

• 2x

t 2 = 1

• p E x

p 1

SLIDE 4

13 cs533d-term1-2005

Discretized potential

Now we have an alternative way to

discretize the equation:

• Approximate the potential energy integral with

a discrete sum

• Take the gradient to get forces

This approach generalizes to all sorts of

forces

Lets do it for multi-dimensional springs

14 cs533d-term1-2005

Multi-dimensional spring potential

Change the L scale in the 1D spring potential to

be the area/volume around the spring

• When we add up, we get an approximation of an

integral over the elastic object

Then get the spring force on i due to j:

Wij = 1

2 E

xi x j Lij 1

• 2

Aij fij = Wij xi = E Aij Lij xi x j Lij 1

• xi x j

xi x j

15 cs533d-term1-2005

Bending

For simulating cloth and thin shells, also

need forces to resist bending

Nontrivial to directly figure out such a force

• n a triangle mesh

Even harder to make sure its roughly

mesh-independent

So lets attack the problem with a potential

energy formulation

16 cs533d-term1-2005

Bending energy

Integrate mean curvature squared over the

surface:

Lets bypass the continuum formulation, and

jump to a discrete approximation of this integral

Split mesh up into regions around “hinges”

(common edges between triangles)

• At each interior edge e, have bending stiffness Be,

curvature estimate e and area of region Ae

Wbend =

1 2 B 2

• Wbend

1 2 Be e 2Ae e:edge

SLIDE 5

17 cs533d-term1-2005

Edge curvature estimate

Look at how the hinge is bent

• Dihedral angle between the incident triangles

Think of fitting a cylinder of radius R parallel to

edge, mean curvature is 1/(2R)

• [side-view of triangle pair: angle between normals is

theta, triangle altitudes are h1 and h2]

For small bend angles (limit as mesh is refined!)

1 2R

• h1 + h2

18 cs533d-term1-2005

Discrete bending energy

As a rough approximation then, using

Ae=|e|(h1+h2)/6:

Treat all terms except theta as constant

(measured from initial mesh), differentiate that w.r.t. positions x Wbend

1 2 Be

e 6Ae

• 2

Ae

e:edge

• =

1 2 Be

e

2

36Ae 2

e:edge

• 19

cs533d-term1-2005

Bending forces

See e.g. Bridson et al. “Simualtion of

clothing…” for reasonable expressions for forces

• Additional simplification: replace theta with a

similarly-behaved trig function that can be directly computed

• Warning: derivation is quite different!