COMS W4167: Physically Based Computer Animation Theme IV, Milestone I - Elastic Simulation Due 22:00:00 PM Tuesday 23th Nov 2021 (EST) Academic Honesty Policy You are permitted and encouraged to discuss your work with other students. Except where explicitly stated otherwise, you may work out equations in writing on paper or a whiteboard. You are encouraged to use Piazza to converse with other students, the TAs, and the instructor. HOWEVER, you may NOT share source code or hard copies of source code. Refrain from activities or the sharing materials that could cause your source code to APPEAR TO BE similar to another student’s enrolled in this or previous years. We will be monitoring source code for individuality. Cheating will be dealt with severely. Source code should be yours and yours only. Do not cheat. For more details, please refer to the full academic honesty policy on the departmental website and on Piazza. 1 Introduction In Theme IV, we’ll tackle the elastic body problem. Simulation of the elastic material is based on formulation of an elastic energy that is quadratic in a strain that describes how the object is deformed. Using different aspects of the deformation as our strain, we can capture different modes of deformation in the energies, and thus obtain the forces that counter-act these modes. In this milestone we’ll examine three elastic forces: spring force, bending force, and the constant-strain-triangle (CST) force. 2 New XML Features This milestone introduces the following new simulation features: 1. The elasticbodyspringforce feature specifies a spring force:
The i1 and i2 attributes specify between which particles this force acts; note that different than the spring force we implemented in the first theme, this force does not require an edge between the two endpoints. The alpha attribute specifies the material’s elastic modulus (equivalent to EA where E is the young’s modulus and A is the area of the cross-section of the spring). The l0 attribute specifies the rest length which means the same as the spring in theme I. 2. The elasticbodybendingforce feature specifies a bending force:
The i1, i2 and i3 attributes specify between which particle this force acts. To represent a bending deformation we need three particles, where the second one is the hinge. The alpha attribute specifies the material’s elastic modulus in resisting the bending mode. The theta0 attribute specifies the rest angle in radian. 3. The elasticbodycstforce feature specifies a CST force:
ebyy="0" ebxy="0" xb1x="0" xb1y="0" xb2x="1" xb2y="0" xb3x="0" xb3y="1"/> 1 Figure 1: The deformation map. The i1, i2 and i3 attributes specify the three vertices of the triangle where this force acts. The youngsmodulus and poisonratio attributes are the material parameters. ebxx, ebxy and ebyy encodes the “resting strain” ¯ that is discussed in the next section. xb1x and xb1y specify the undeformed position of vertex 1, and xb2x, xb2y, xb3x and xb3y are similar. 3 Kinematics 3.1 The Deformation Map and Deformation Gradient The kinematics of a deforming material are represented by a deformation mapping Φ : R2 → R2 taking every ‘undeformed’ material point x on the plane to its corresponding deformed position Φ(x). In general, Φ is an arbitrary nonlinear map taking 2 coordinates denoting a position to a new set of 2 coordinates denoting some other position. Now we’d like to know how much this deformation is stretching, compressing or shearing the material. Since a generic deformation is an arbitrary nonlinear map, we only look at each point locally (by examining its infinitesimal neighborhood). Consider a point x¯ in the undeformed space, and a nearby point x¯ + δx¯. The deformed position are x and x + δx, and by definition of Φ we have: x = Φ(x¯) x + δx = Φ(x¯ + δx¯) = Φ(x¯) +∇Φ(x¯)δx¯ +O(‖δx¯‖2) δx = ∇Φ(x¯)δx¯ +O(‖δx¯‖2) ∇Φ(x¯) is the Jacobian of the transformation Φ evaluated at x¯, often called the Deformation Gradient. The high order terms in the end can be omitted because we assume δx¯ is infinitesimal. 3.2 The Strain Now we are ready to evaluate how much lengths change after such a transformation. Under the assumption that δx¯ is infinitesimal, the line from x¯ to x¯+δx¯ in the undeformed space remains a line after the deformation, 2 going from x to x+ δx. Therefore the lengths of this piece of material before and after deformation are ‖δx¯‖ and ‖δx‖ respectively. For the ease of computation we look at the change in squared length (so as to avoid taking square roots): δxT δx− δx¯T δx¯ = (∇Φ(x¯)δx¯)T (∇Φ(x¯)δx¯)− δx¯T δx¯ = δx¯T∇Φ(x¯)T∇Φ(x¯)δx¯− δx¯T δx¯ = δx¯T [∇Φ(x¯)T∇Φ(x¯)− I]δx¯ The reason we can use change in squared length instead of change in length itself as a measure of stretching is that, to the first order, they are equivalent as long as the deformation is small. In order to make this change in squared length a scale-invariant representation of relative stretching/compression, we extract out the term in the middle, [∇Φ(x¯)T∇Φ(x¯)− I], and call it the strain, or . This strain is a second-order tensor, unlike the 1D case where a single scalar fully describes how much a spring is stretched as shown in class two weeks ago. In 2D materials can be stretched differently in different directions, so a second-order tensor is needed. The physical meaning of this tensor is that when you hit it on both sides with a direction δx¯, it returns a change in squared length in that direction. In its matrix representation, the eigenvectors of this matrix are the two directions that are stretched or compressed the most, and the corresponding eigenvalues tell you how much the stretching/compression is in that direction. Another interesting thing to note is that is rotation-invariant, even though ∇Φ is not. As a measure of deformation, the strain should not depend on any rigid motion such as translation and rotation. It is obvious that ∇Φ, which is the result of a spatial differential operator, stays the same under translation, but what about rotation? To test this, we apply a constant rotation Q to any deformation Φ to obtain a new deformation Φ˜, and see how responds: Φ˜ = QΦ ∇Φ˜ = Q∇Φ ˜ = ∇Φ˜T∇Φ˜− I = (Q∇Φ)T (Q∇Φ)− I = ∇ΦTQTQ∇Φ− I = ∇ΦT∇Φ− I = Thus we’ve convinced ourselves that our strain is invariant under any rigid motions, therefore it qualifies as a measure of deformation. 3.3 Example Deformation Let’s look at an example with real numbers. Consideration an affine deformation of a rectangular sheet of rubber as in Figure 2. This rectangle is stretched along the green axis by 80%, and compressed along the red axis by 40%. This means our deformation map Φ is multiplication by a matrix with two eigenvalues: 1.8 corresponding to the eigenvector (3, 1) (green axis), and 0.6 corresponding to the eigenvector (-1, 3) (red axis). Therefore Φ looks like this: Φ(x) = 1√ 10 ( 3 −1 1 3 )( 1.8 0 0 0.6 ) 1√ 10 ( 3 1 −1 3 ) x = ( 1.68 0.36 0.36 0.72 ) x 3 Figure 2: An example deformation of square, with stretching in one direction and compression in the other. The two arrows denote the two eigenvectors of the strain tensor with corresponding eigenvalues (λ1 and λ2) marked on the arrows. We can compute the strain tensor from Φ: = ∇Φ(x¯)T∇Φ(x¯)− I = ( 1.68 0.36 0.36 0.72 )T ( 1.68 0.36 0.36 0.72 ) − I = ( 1.952 0.864 0.864 0.352 ) It is easy to verify that this 2x2 matrix has eigenvalues 2.24 and -0.64, corresponding to a 80% stretching (1.82 − 1 = 2.24) and a 40% compression (0.62 − 1 = −0.64). 4 Physics The strain tensor gives us a complete description of the deformation of the material, but it doesn’t tell us how the material will respond. That is out of the scope of kinematics, and we have to put in some ingredients from physics. Here in this section we’ll attempt to write down the energy that arises from the deformation, which governs how the simulation will take place in time. For this purpose, we can proceed in either a mechanics-oriented way, or a mathematics-oriented way. 4.1 The Mathematical Way Imagine we have somehow obtained the energy function w of strain , we can take its Taylor expansion: w() = w(0) + ∂w(0) ∂ + 1 2 ∂2w(0) ∂2 2 +O(‖‖3) w(0) is a constant independent of and can be conveniently assume to be 0. The linear term is also 0 because the energy should be minimized at zero strain (the definition of “resting shape”). If we adopt the small strain assumption, the higher order terms can be omitted since the energy is dominated by the quadratic term. This is basically saying, no matter how complicated the material laws are, under the small strain assumption we can always write down the energy as a quadratic function of the strain. 4.2 The Mechanical Way Physically, the energy is equal to the product between the stress and the strain: 4 w() = 1 2 2∑ i=1 2∑ j=1 ijσij and the stress σ is defined as the elastic moduli times the strain: σkl = 2∑ i=1 2∑ j=1 cijklij therefore w() = 1 2 2∑ i=1 2∑ j=1 2∑ k=1 2∑ l=1 cijklijkl Here we obtain an energy that is quadratic in strain again. c is a fourth-order tensor that encodes all the material properties such as stiffness upon stretching/compression in each direction and shearing etc. Although a generic fourth-order tensor requires a 16 scalars to represent, symmetry in both the strain tensor and the stress tensor reduces this number to 9. If we further assumes that the material is isotropic, then we only need two parameters to fully describe this material: Young’s modulus E and Poisson ratio ν: σxxσyy σxy = E 1− ν2 1 ν 0ν 1 0 0 0 1− ν 2 xxyy 2xy 4.3 Total Energy The w() defined above gives the energy density at a point in the material. Since is a function of undeformed space position x¯, so is w. To obtain the total elastic energy that is stored in an elastic body, we integrate this energy density over the entire body: W = ∫ Ω w(x¯)dΩ For this milestone, we simply assume that the strain is constant for any element (discussed in the next section), therefore the energy is simplified to W = w(x¯)A where A is the 2D element’s area. 4.4 Discretization All the analysis up to now has been performed on a continuum, in the form of differential equations. Continua have an infinite number of degrees of freedom and cannot be manipulated directly in computers - we need to discretize. The discretization here is in the spatial dimensions, but it’s in the same spirit with the temporal dis- cretization we’ve implemented in the first Theme with integrators. By spatial discretization we represent the body of interest with a set of nodes, connected by edges, polygon faces or polyhedron volumes, depending on the dimension of the body. These building blocks are called elements. There are a fixed number of nodes in each element, and due to the limited amount of information these nodes can provide, we have to assume 5 certain properties on the element itself. For example, a triangle element with only one node in each corner (one of the simplest and most commonly used element types) simply does not have any information about how the deformation mapping changes over the triangle to the second or higher order, so we have to assume that the deformation Φ is linear, and when we need information about a point inside the triangle, we linearly interpolate between the three nodes. Linear deformation map implies constant deformation gradient, and as a result this element type is called constant strain triangle (discussed later in more details). Of course the interpolation won’t be able to faithfully reproduce the true property values at an arbitrary point in the element, so a discretization error is introduced; but since our hardware technology isn’t advanced enough to handle infinite degrees of freedom, this error is inevitable. As long as the discretization error approaches zero as we refine our mesh further and further, we say the discretization converges and thus is valid. For our milestone, we do not consider more complicated (but potentially more accurate) elements. 5 Required Features 5.1 Elastic Body Spring Force The spring force is similar to the spring force we implemented in Theme 1, with the only difference being that we use a shape-independent material property as the stiffness. Following the concept from above, we define the strain as: = ‖xi − xj‖ − l0 l0 The energy density is: w = 1 2 α2 = α 2l20 (‖xi − xj‖ − l0)2 Assuming that the strain is constant along the spring, the total energy is W = wl0 = α 2l0 (‖xi − xj‖ − l0)2 Taking the gradient of this energy gives the force. Please implement the computation for energy, force, and hessian in addEnergyToTotal(), addGradEToTotal() and addHessEToTotal() in ElasticBodySpringForce.cpp. Note that implementation of the energy is optional; it won’t be graded. Implementing the energy will help implementing the force by providing more intuitions and more debugging information. 5.2 Elastic Body Bending Force When a thin rod is bent, why does it try to restore its initial straight configuration? It is because the bending introduces differential compression/stretching along the cross section of the rod, and thus increases the elastic energy. If we explicitly model the cross section of the rod, we will have an array of axial springs side-by-side, connected by short and highly stiff radial-directed springs at endpoints. When bent springs on the inner layer are compressed and springs on the outer layer are stretched, therefore creating the restoring force. However this modeling scheme has two serious disadvantages: it is too stiff and has bad aspect ratios, both of which lead to numerical difficulties. In practice we model thin rods as a single polyline, augmented with a bending force to resist bending. Since the bending energy arises from axial compression and stretching of the material, it is naturally quadratic in the curvature κ of the material centerline: 6 Figure 3: The bending at the common node between two edges w = α 2 κ2 where α captures all other constants, including the material’s elasticity modulus and cross section area. As our discretization refines, curvature κ can be approximated by θ ‖e¯ij‖+ ‖e¯jk‖ , where θ is the angle by which the next edge (eij) “deviates” from the direction of the previous edge (ejk): θ = atan2(eij × ejk, eij · ejk) where function atan2(y, x), a C library function, returns the direction angle of vector 〈x, y〉. Therefore the total bending energy is: W = w(‖e¯ij‖+ ‖e¯jk‖) = α 2 κ2(‖e¯ij‖+ ‖e¯jk‖) ≈ α 2 ( θ ‖e¯ij‖+ ‖e¯jk‖ ) 2(‖e¯ij‖+ ‖e¯jk‖) = α 2(‖e¯ij‖+ ‖e¯jk‖)θ 2 Allowing non-straight (non-zero curvature or non-zero bending angle) resting shape, we define the energy as W = α 2(‖e¯ij‖+ ‖e¯jk‖) (θ − θ0) 2 Please calculate the gradient and hessian, and implement addEnergyToTotal(),anddGradEToTotal() and addHessEToTotal() in ElasticBodyBendingForce.cpp. Note that implementation of the energy is optional; it won’t be graded. Implementing the energy will help implementing the force by providing more intuitions and more debugging information. The following formula may be useful: 7 Figure 4: A triangle undergoing a deformation Φ ∂arctan(x) ∂x = 1 1 + x2 ( ∂a× b ∂a )T = ( 0 1 −1 0 ) b 5.3 Elastic Body Constant Strain Triangle Force Both the spring force and the bending force essentially assumes an infinitely thin rod-like structure, where the strain is fully captured by a single scalar (either the edge length change or the bending angle). For 2D elastic bodies we need to use the full-fledged treatment in the previous sections. Now consider the triangle involving particle i, j and k: We assume that the strain is constant over the interior of the triangle (hence the name “constant strain triangle”), and so the deformation map Φ is affine (without loss of generality, from here on we assume i = 1, j = 2 and k = 3): x = ∇Φ(x¯− x¯1) + x1 ∇Φ = ( x2 − x1 x3 − x1 y2 − y1 y3 − y1 )( x¯2 − x¯1 x¯3 − x¯1 y¯2 − y¯1 y¯3 − y¯1 )−1 where xi denotes the x component of node xi, yi denotes its y component, and the overhead bar denotes the undeformed configuration as before. Denoting the components in ∇Φ as ∇Φ = ( a c b d ) , we define the strain: 8 = ( xx xy xy yy ) = ∇ΦT∇Φ− I = ( a2 + b2 − 1 ac+ bd ac+ bd c2 + d2 − 1 ) and the energy: W = A¯ 2 σxxσyy σxy T xxyy 2xy = A¯ 2 ( E 1− ν2 1 ν 0ν 1 0 0 0 1− ν 2 xxyy 2xy )T xxyy 2xy = EA¯ 2(1− ν2) [ 2 xx + 2νxxyy + 2 yy + 2(1− ν)2xy] where A¯ is the undeformed triangle area: A¯ = 1 2 ‖(x¯3 − x¯1)× (x¯2 − x¯1)‖ In order to find the force, we take the gradient of the energy above: −F = ∇xW = EA¯ 1− ν2 [(xx + νyy)∇xxx + (yy + νxx)∇xyy + 2(1− ν)xy∇xxy] ∇xxx = 2a∇xa+ 2b∇xb ∇xyy = 2c∇xc+ 2d∇xd ∇xxy = a∇xc+ c∇xa+ b∇xd+ d∇xb Now letting ( x¯2x − x¯1x x¯3x − x¯1x x¯2y − x¯1y x¯3y − x¯1y )−1 = ( e g f h ) we have ∇Φ = ( (−e− f)x1 + ex2 + fx3 (−g − h)x1 + gx2 + hx3 (−e− f)y1 + ey2 + fy3 (−g − h)y1 + gy2 + hy3 ) ∇a = −e− f 0 e 0 f 0 ,∇b = 0 −e− f 0 e 0 f ,∇c = −g − h 0 g 0 h 0 ,∇d = 0 −g − h 0 g 0 h , 9 With these derivation, you should be able to implement the elastic energy, force, and hessian of a constant strain triangle. Please do so in addEnergyToTotal(), addGradEToTotal(), and addHessEToTotal() in ElasticBodyCSTForce.cpp. Note that implementation of the energy is optional; it won’t be graded. Implementing the energy will help implementing the force by providing more intuitions and more debugging information. Debugging help If you have problems implementing the elasticity, try debugging a few simple cases where you know the exact solution. Then test the output of your algorithm for the deformation of each point to the exact solution that you know to be true. An easy example to use for this is the triangle with vertices (0, 0), (0, 1) and (1, 0) (ref. the first image below). Another simple test case is the square with vertices (0, 0), (0, 1), (0, 1) and (1, 1) (ref. the second image below) A triangle with vertices (0, 0), (0, 1) and (1, 0) A square with vertices (0, 0), (0, 1), (0, 1) and (1, 1) 10 6 Creative Scene As part of your final submission for this milestone, please include a scene of your design that best shows off your program. Based on the quality of your scene, you will have the opportunity to earn up to 10% extra credit. Your scene will be judged by a secret committee of top scientists using the highly refined criteria of: 1. How well the scene shows off this milestone’s ‘magic ingredients’ (a la Iron Chef ). 2. Aesthetic considerations. The more beautiful, the better. 3. Originality. Top examples will be posted to Piazza. To submit this scene, place the XML file in the assets/Creative/ directory of your submission, and submit a video of the scene to Courseworks. To export PNGs of your scene, run path/to/FOSSSim -s path/to/creativescene.xml -m path/to/output/dir/ This will output a PNG file for each frame of the simulation. The program will output black frames if you pass it -d 0 so leave that flag out. Then, create a video from your PNGs using ffmpeg ffmpeg -r 24 -f image2 -i pngs/frame%05d.png -vcodec libx264 -pix_fmt yuv420p out.mp4 If you are having issue with saving using svgs, please refer to Piazza post 6.1 Creative scene submission Please download and submit the mp4 file to Courseworks. Remember to generate videos of your creative scene before submitting, as the Codio box will become read-only after you submit. The due date for creative scene will be on 12:00pm (noon) Nov. 24 Please refer to announcement on Courseworks 11 欢迎咨询51作业君