Notes Modern FEM � Galerkin framework (the most common) � Assignment 2 is up � Find vector space of functions that solution (e.g. X(p)) lives in • E.g. bounded weak 1st derivative: H 1 � Say the PDE is L[X]=0 everywhere (“strong”) � The “weak” statement is � Y(p)L[X(p)]dp=0 for every Y in vector space � Issue: L might involve second derivatives • E.g. one for strain, then one for div sigma • So L, and the strong form, difficult to define for H 1 � Integration by parts saves the day cs533d-term1-2005 1 cs533d-term1-2005 2 Weak Momentum Equation Making it finite � The Galerkin FEM just takes the weak equation, � Ignore time derivatives - treat acceleration and restricts the vector space to a finite- as an independent quantity dimensional one • We discretize space first, then use • E.g. Continuous piecewise linear - constant gradient “method of lines”: plug in any time integrator over each triangle in mesh, just like we used for Finite Volume Method [ ] = � ˙ X � f body � � � � ˙ L X � This means instead of infinitely many test functions Y to consider, we only need to check a ( ) � � [ ] = Y � ˙ X � f body � � � � ˙ finite basis Y L X � � � The method is defined by the basis � � � = Y � ˙ ˙ � � Y � � � • Very general: plug in whatever you want - X Yf body polynomials, splines, wavelets, RBF � s, … � � � � � � = Y � ˙ ˙ � + � � Y X Yf body � � � cs533d-term1-2005 3 cs533d-term1-2005 4
Linear Triangle Elements The equations � Simplest choice � Take basis { � i } where � � � � � j � ˙ i � i � � j f body + � � � j = 0 x ˙ � i (p)=1 at p i and 0 at all the other p j � s i • It � s a “hat” function � � � � Then X(p)= � i x i � i (p) is the continuous piecewise linear � � � � �� j � i ˙ = � j f body � � � � j ˙ x function that interpolates particle positions i i � Similarly interpolate velocity and acceleration � � � � Plug this choice of X and an arbitrary Y= � j (for any j) into •Note that � j is zero on all but the triangles the weak form of the equation � Get a system of equations (3 eq. for each j) surrounding j, so integrals simplify •Also: naturally split integration into separate triangles cs533d-term1-2005 5 cs533d-term1-2005 6 Change in momentum term Body force term � m ij = �� i � j � Usually just gravity: f body = � g � Let � � Rather than do the integral with density all over m ji ˙ x ˙ � Then the first term is just again, use the fact that � I sum to 1 i i • They form a “partition of unity” M ˙ ˙ x � Let M=[m ij ]: then first term is • They represent constant functions exactly - just about � M is called the mass matrix necessary for convergence • Obviously symmetric (actually SPD) � Then body force term is gM1 • Not diagonal! � More specifically, can just add g to the � Note that once we have the forces (the other accelerations; don � t bother with integrals or integrals), we need to invert M to get accelerations mass matrix at all cs533d-term1-2005 7 cs533d-term1-2005 8
Stress term The algorithm � Calculate constant strain and strain rate (so � Loop over triangles constant stress) for each triangle separately • Loop over corners � Note �� j is constant too • Compute integral terms � So just take ��� j times triangle area � only need to compute M once though - it � s constant • End up with row of M and a “force” � [derive what �� j is] � Solve Ma=f � Magic: exact same as FVM! • In fact, proof of convergence of FVM is often (in other � Plug this a into time integration scheme settings too) proved by showing it � s equivalent or close to some kind of FEM cs533d-term1-2005 9 cs533d-term1-2005 10 Lumped Mass Consistent vs. Lumped � Inverting mass matrix unsatisfactory � Original mass matrix called “consistent” • For particles and FVM, each particle had a mass, so � Turns out its strongly diagonal dominant (fairly we just did a division easy to solve) • Here mass is spread out, need to do a big linear solve - even for explicit time stepping � Multiplying by mass matrix = smoothing � Idea of lumping: replace M with the “lumped � Inverting mass matrix = sharpening mass matrix” � Rule of thumb: • A diagonal matrix with the same row sums • Implicit time stepping - use consistent M • Inverting diagonal matrix is just divisions - so diagonal (counteract over-smoothing, solving system anyways) entries of lumped mass matrix are the particle • Explicit time stepping - use lumped M masses • Equivalent to FVM with centroid-based volumes (avoid solving systems, don � t need extra sharpening) cs533d-term1-2005 11 cs533d-term1-2005 12
Locking Locking and linear elements � Simple linear basis actually has a major problem: � Look at nearly incompressible materials locking � Can a linear triangle mesh deform incompressibly? • But graphics people still use them all the time… • [derive problem] � Notion of numerical stiffness � Then linear elements will resist far too much: numerical • Instead of thinking of numerical method as just getting stiffness much too high an approximate solution to a real problem, � Numerical material “locks” • Think of numerical method as exactly solving a � FEM isn � t really a black box! problem that � s nearby � Solutions: • For elasticity, we � re exactly solving the equations for a • Don � t do incompressibility material with slightly different (and not quite • Use other sorts of elements (quads, higher order) homogeneous/isotropic) stiffness • … � Locking comes up when numerical stiffness is MUCH higher than real stiffness cs533d-term1-2005 13 cs533d-term1-2005 14 Quadrature Accuracy � Formulas for linear triangle elements and constant � At least for SPD linear problems (e.g. density simple to work out linear elasticity) FEM selects function from � Formulas for subdivision surfaces (or high-order finite space that is “closest” to solution polynomials, or splines, or wavelets…) and varying • Measured in a least-squares, energy-norm density are NASTY � Instead use “quadrature” sense • I.e. numerical approximation to integrals � Thus it � s all about how well you can � Generalizations of midpoint rule approximate functions with the finite space • E.g. Gaussian quadrature (for intervals, triangles, tets) or tensor you chose products (for quads, hexes) � Make sure to match order of accuracy [or not] • Linear or bilinear elements: O(h 2 ) • Higher order polynomials, splines, etc.: better cs533d-term1-2005 15 cs533d-term1-2005 16
Hyper-elasticity Variational Derivatives � Another common way to look at elasticity � Force is the negative gradient of potential • Useful for handling weird nonlinear compressibility laws, for • Just like gravity reduced dimension models, and more � What does this mean for a continuum? � Instead of defining stress, define an elastic potential • W=W( � X/ � p), how do you do -d/dX? � � W � X � p + �� Y energy � [ ] = � Variational derivative: W total X + � Y � � • Strain energy density W=W(A) � � p � • W=0 for no deformation, W>0 for deformation � � • So variational derivative is W � X � + �� W � Y � � • Total potential energy is integral of W over object � - � • � W/ � A � p � A � p � � � This is called hyper-elasticity or Green elasticity • And f= � • � W/ � A � W � Y � • Then stress is � W/ � A = W total + � � For most (the ones that make sense) � A � p stress-strain relationships can define W � Easy way to do reduced Y � � � W � • E.g. linear relationship: W= � : � =trace( � T � ) dimensional objects = W total � � � A (cloth etc.) cs533d-term1-2005 17 cs533d-term1-2005 18 Numerics � Simpler approach: find discrete W total as a sum of W � s for each element • Evaluate just like FEM, or any way you want � Take gradient w.r.t. positions {x i } • Ends up being a Galerkin method � Also note that an implicit method might need Jacobian = negative Hessian of energy • Must be symmetric, and at least near stable configurations must be negative definite cs533d-term1-2005 19
Recommend
More recommend
