Learning Elastic Constitutive Material and Damping Models
Bin Wang, Yuanmin Deng, Paul Kry, Uri Ascher, Hui Huang, Baoquan Chen
11 Learning Elastic Constitutive Material andDamping Models
Bin Wang Paul Kry Yuanmin Deng Uri Ascher Hui Huang Baoquan Chen
Abstract —The fidelity of a deformation simulation is highly dependent upon the underlying constitutive material model. Commonlyused linear and nonlinear constitutive material models contain many simplifications and only cover a tiny part of possible materialbehavior. In this work we propose a framework for learning customized models of deformable materials from sparse example surfacetrajectories. The key idea is to iteratively improve a correction to a nominal model of the elastic and damping properties of the object,which allows new forward simulations with the learned correction to more accurately predict the behavior of a given soft object. Thechallenge is that such data is sparse as it is typically available only on part of the surface. Sparse reduced space-time optimizationidentifies gentle control forces with which we extract necessary annotated data for model inference and to finally encapsulate thematerial correction into a compact parametric form. We demonstrate our method with a set of synthetic examples, as well as with datacaptured from real world homogeneous elastic objects.
Index Terms —Computing methodologies, physical simulation, elastic and damping model, dynamics (cid:70)
NTRODUCTION T HE simulation of deformable objects is ubiquitous incomputer graphics and robotics research due to thelarge number of varied applications, including animation,movie making, medical treatment and manufacturing. Toobtain accurate deformation simulation, finite element mod-elling and continuum elasticity laws are typically employed,where a constitutive material model that can cover the rangeof material behavior is selected. Subsequently, the materialparameters are needed to be carefully tuned in order to fitempirical data. Data driven based methods have recently ex-hibited great potential in this direction. Advanced scanningtechnologies can be used to faithfully capture a deformationbehavior under external force, and in turn the data can beused to estimate the parameters of the mathematical model.However, currently there is no standard solution forchoosing appropriate constitutive models, especially forlarge deformations, real-world material modelling, fine levelsimulation of heterogeneous models, and artist designedcartoon physics. The situation for modelling damping iseven worse. Indeed, there is no agreement about dampingmodel in the mechanical engineering literature which candescribe various damping effects in a unified way. Manyuse the Rayleigh model, but it can be inadequate for visualpurposes [1]. Based on these observations, in this paper • Bin Wang is with Advanced Innovation Center for Future Visual Enter-tainment of Beijing Film Academy. E-mail: [email protected] • Paul Kry is with School of Computer Science, McGill University. E-mail:[email protected] • Yuanmin Deng is with School of Computer Science and Technology,Shandong University. E-mail: [email protected] • Uri Ascher is with Department of Computer Science, University of BritishColumbia. E-mail: [email protected] • Hui Huang is with College of Computer Science and Software Engineer-ing, Shenzhen University. E-mail: [email protected] • Baoquan Chen is with Advanced Innovation Center for Future VisualEntertainment of Beijing Film Academy and Peking University. E-mail:[email protected] • Hui Huang and Baoquan Chen are the corresponding authors of this paper. we propose a more general material inference frameworkobtained by directly learning constitutive laws rather thanfitting parameters. Similar ideas have also been used by [2]and [1], which focus on providing a forward modeling toolfor material design. Learning a constitutive model is chal-lenging because: (1) the model should be able to encapsulateall the variations of material properties in a generic way; (2)there may be no obvious source of training data.The main contribution of this paper are:(1) The constitutive material model is designed as acombination of empirical baseline model and a paramet-ric representation, in which all the variations of materialproperties are encapsulated. This model maintains a goodbalance between stability and flexibility.(2) We propose a progressive inverse learning frame-work, which is capable to learn a complex constitutivematerial from sparse motion trajectories. The highlight inthis framework is leveraging a space-time optimizationtechnique to generate annotated data for constitutive modelfitting.(3) By incorporating a probabilistic correspondence rep-resentation, a differentiable shape matching metric is devel-oped. Significantly, this allows our system to work faithfullyon real-world captured data.Fig. 1 shows a preview of our approach and results. Wedemonstrate the performance of our approach on severalproblems, including synthetic examples, coarsening appli-cations, and captured data. The variety of results describedin Section 7 leads us to conclude that we have achievedan important step towards the data based constructionand understanding of nonlinear constitutive elastodynamicforce models. a r X i v : . [ c s . G R ] S e p Fig. 1: Our parametric material model learns a correction to a nominal material model from kinematic data alone, allowingus to accurately capture the nonlinearity of different constitutive material models. Left: classical nonlinear constitutivematerial. Middle: user designed elasticity and damping. Right: real world material.
ELATED WORK
With recent improvements in sensing technologies, the data-driven approach of modeling and reconstructing deforma-tion parameters from real world measurements has offeredgreat potential for computer graphics applications, such asfabrics, soft objects, and human organs and faces [3], [4],[5], [6], [7], [8]. Bickel et al. [8] fit material parameters withan incremental loading strategy to better approximate non-linear strain-stress relationships. Wang et al. [6] proposeda piecewise linear elastic model to reproduce nonlinear,anisotropic stretching and bending of cloth. Miguel et al.[7] directly optimized nonlinear stress-strain curves basedon measurements. Then Miguel et al. [9] estimated internalfriction. Further, Miguel et al. [10] developed a method formodeling example based materials with energy functionsfor both cloth and elastic solids.A common limitation with previous methods is thatthey require a dense force displacement field. While Bhatet al. [11] avoided the need for force capture by using videotracking of cloth, they still assumed a trivial cloth referenceshape. Yang et al. [12] presented a learning-based algorithmto recover material properties of cloth from videos, usingtraining data sets generated by physics-based simulators,and their focus was on material type estimation due toinconsistency between real and synthetic data and sparsematerial space sampling. Davis et al. [13] estimated materialand damping properties through extracting small vibrationmode from high-speed and regular frame-rate video. Wanget al. [14] estimated linear elastic material parameters frompartially observed surface trajectories of an object’s passivedynamics. Our work has a similar setting, but focuses oncorrecting the errors that arise from assuming simple elasticand damping force models.
Another popular trend in computer graphics is to directlyfit parametric functions as a material description. Xu et al.[2] provided a method to design isotropic and anisotropic(orthotropic) nonlinear solid elastic materials using a piece-wise spline interface, which can provide local control ondeformation behavior. A hyperelasticity model based on en-ergy addends proposed by [10] allows modelling of variousnonlinear elasticity effects in a separable manner. Instead of explicitly modelling the stress-strain relationship, Martin etal. [15] and Schumacher et al. [16] promoted an art-directedapproach to solid simulation, which constructs a manifoldof prefered deformations by examples and guides the objecttowards it. Coarsening techniques have proved useful forcomputational design for fabrication [17], [18], [19], whereone must deal with the problem of creating equivalentphysics based models. Kharevych et al. [20] took an en-ergy based approach to coarsening composite elastic objectsthrough the use of global harmonic displacements. Nesmeet al. [21] created nonlinear shape functions and projectedfine-level mass, stiffness, and damping matrices to producecoarse composite elements, while Torres et al. [22] intro-duced an improved element based coarsening method thatdeals with co-rotation. Compared with elasticity modelling,few publications have focused on the design of dampingmodel, the work of [1] being the first to present a designmethod for anisotropic and/or nonlinear damping effects.Banderas et al. [23] focuses on dissipation potentials baseddamping modelling, similarly to us using strain rate tocontrol damping. Targeting cloth hysteresis effects, Miguelet al. [9] proposed a internal friction model based on anaugmented reparameterization of Dahl’s model.In our work, inspired by [2] we use principal stretchesto formulate nonlinear elastic and damping force correctionsdue to its complexity and flexibility. However, instead of us-ing a spline interface, we adopt RBF and NN representationsto parameterize the material correction.
NN analysis has had a liberating effect on materials science,by enabling the study of incredibly diverse phenomenawhich are not as yet accessible to physical modelling [24],[25]. Jung and Ghaboussi [26] modeled rate-dependent ma-terials with NN, giving results both for a synthetic exampleand for data from a pre-stressed concrete beam. Stefanos andGyan [27] used the length of strain trajectory traced by a ma-terial point, also called intrinsic time, as an additional inputparameter in training. This is essential for situations of cyclicand transient loading. To overcome the challenge of traininga convolutional NN with a small dataset, Liang et al. [28]employed a training strategy that combines three key ideas:unsupervised deep learning to determine the filter parame-ters of a convolution layer (generally using encoder-decoder based unsupervised learning strategies), supervised learn-ing to determine the parameters in the classifier or regressorlayer, and data augmentation to generate more trainingdata. In the computer graphics community, deep learningtechnology for deformation modelling has gradually gainedmore attention. The DeepWarp technique of [29] attemptsto learn a mapping from a linear elasticity simulation toits nonlinear counterpart. Fulton et al. [30] perform timeintegration of the elastodynamic system in a learned non-linear reduced latent-space, which is represented using anartificial NN. In comparison, the algorithm proposed in ourpaper puts emphasis on how to accurately approximate theunderlying constitutive model for FEM simulation. We willdemonstrate that traditional computer graphics algorithmscan be adapted to address an interesting and timely probleminvolving data labeling and augmentation.
VERVIEW
The core of our approach is to learn accurate elastic anddamping properties of an object through trajectory fitting. Abaseline constitutive elastodynamic model, called nominalmodel in Fig. 2, is assumed to be given The inferred baselinemodel correction is encapsulated in RBF or Neural Networkparametric representation.The input of the system is a set of incomplete surfacetrajectories of an object moving dynamically, unforced, inresponse to an initial perturbation. The main loop alternatesbetween solving a sparse reduced space-time optimizationproblem, and fitting a correction function (see Fig. 2). Thesparse reduced space-time optimization softly constrainsnodes of the tetrahedral mesh to follow a sparse collec-tion of incomplete surface trajectories in addition to thediscretized physical equations of motion. Through injectinggentle forces we can keep the system trajectories close todesired input trajectories and still obey Newton’s physicallaws. The key idea is that this gentle control force identifieswhat is currently missing due to material model inaccuracyand should be compensated using the correction model. Ad-ditionally, we solve an overdetermined problem to identifythe missing corrective stress on each tetrahedron from thevertex control forces taking all the frames into consideration.Then, the best fit correction model is distilled from thestrain, strain rate, and stress data. In the following iterations,the current correction of the nominal material model isadded to both the forward simulation and the space-timeoptimization. As the iterations progress, the correction isgradually refined to provide better accuracy, which in turnimproves how the forward simulation matches the exampletrajectory. During each optimization iteration, the correctionis encapsulated in an RBF parametric representation, and itis only after the final iteration that we train a compact NNbased representation with data sampled from the final RBFfunction.At the beginning of each iteration, we complete a for-ward simulation with the updated material model to eval-uate how well this trajectory matches the example surfacetrajectory to determine convergence. We can similarly mon-itor the magnitude of gentle control forces identified bythe space-time optimization at each loop, and continue toiterate as long as we see improvement, even if the forward simulation error alone does not reveal that progress is beingmade. We also use this regular forward simulation trajectoryto extract a reduced basis for the optimization. Moreover,We compute an initial seed trajectory for this simulation us-ing a forward simulation, which includes constraints on theimmobilized parts of the mesh, as well as constraints to en-sure that the trajectories of the sparse surface points followtheir known positions. The constrained forward simulationgives a good starting point for the sparse reduced space-time optimization, which quickly converges to a solutionthat identifies a plausible trajectory for unobserved nodesand the corresponding gentle control forces.
ATERIAL M ODEL
Even though there are plenty of empirical hyperelastic con-stitutive material models such as the nonlinear St. Venant-Kirchhoff, neo-Hookean, Ogden or Mooney-Rivlin materi-als, they do not account for all deformation phenomena thatmay arise. Choosing a correct model to fit measurement datais already a difficult task. Our data-driven approach allowsus to learn a parametric material model that can encapsulatea wide range of elastic and damping properties in a compactand unified correction function.
Our deformable models are constructed using linear shapefunctions. In order to handle large deformations of softobjects, the nominal material is described in terms of thewidely adopted co-rotated linear FEM, formulated usingprincipal stretches [2]. The ensuing computation of theelement stresses and vertex forces is straightforward.The deformation gradient F for each tetrahedron isdiagonalized by SVD, F = U ˆ F V T , and the first Piola-Kirchoff stress is computed with the principal stretches, ˆ P ( ˆ F ) = 2 µ ( ˆ F − I ) + λtr ( ˆ F − I ) I where µ and λ areLam´e parameters. The diagonal stress is then transformedback to the world frame, P = U ˆ P ( ˆ F ) V T . An element’scontribution to its vertex forces is
P B m , where B m is theinverse material space shape matrix (see [31]). Summing thecontribution of all elements, we can build a large sparsematrix B , which combines the entries in U , V , and B m , andcan be multiplied by the block vector of all element diagonalstresses ˆ p to give a block vector of all vertex forces f , that is, B ˆ p = f .We also include Rayleigh damping in our nominalmodel, with forces computed by ( α M + α K )˙ x , where ˙ x are the FEM vertex velocities, M is the lumped massmatrix, and K is the stiffness matrix assembled from per-element stiffness matrices. The nominal model parametersare assigned manually or computed with the method of [14]. Our parametric material correction model ultimately com-putes a correction to the Piola stress to compensate forthe elasticity and damping simultaneously. Similar to [32],the stress correction is computed in a rotated frame fromdiagonalized strain and strain-rate as ∆ ˆ P = N ( ˆ F , ˆ˙ F ) , andthen rotated back to current world frame by the same U and V used to diagonalize F . Here, ˆ˙ F is defined as the 3 Fig. 2: Schematic overview of our data-driven parametric material learning framework which iteratively learns a correctionto a nominal material model that allows us to accurately reproduce the captured trajectory, even when the nominal modeldiffers significantly.diagonal components of U ˙ F V T , where ˙ F is the deformationgradient velocity. Using principal stretches and the diagonaldeformation gradient velocity reduces the complexity of ourfunction approximation problem, but still permits complexstress corrections and strain dependent damping.As discussed later, the 6-input 3-output function approx-imation problem allows us to use an RBF with a moderatenumber of basis functions, and select a network architecturewith a moderate number of hidden nodes. The correctedstresses in the world frame are P n = U ( ˆ P ( ˆ F ) + N ( ˆ F , ˆ˙ F )) V T . (1)For implicit integration, the gradient of N can be easilycomputed from the RBF, or from the NN function with auto-matic differentiation; the gradient of stress P with respect tothe deformation gradient F can be computed by the productrule and a careful evaluation of the different terms [2]. ORCE C ORRECTION E STIMATION
The desired surface trajectory provides a rich source of in-formation about the dynamics of the object. The purpose ofusing space-time optimization is to compute a set of gentlecontrol forces that will drive the simulation follows captureddata. Concequently, nessesary information to correct ourcurrently estimated parametric material model such thatthe simulation follows captured data can be distilled fromspace-time optimization result .
Many variations of the space-time constraints approach of[33] have been proposed. To deal with the large numberof degrees of freedom in our deformable models, we usereduction and sparse constraints taking inspiration from[34] and [35].We compute a reduce basis Φ using a data driven basedmethod Proper Orthogonal Decomposition(POD), whichperform principal component analysis (PCA) on forwardsimulation trajectory with the provided initial conditionsand our current approximation of the material model. Forthe analysis we use a short portion of the forward simula-tion sequence which targets the most interesting dynamics,and perform PCA on only a fraction of the frames, in turnkeeping only a fraction of the vectors for the basis. We solve the space time optimization as an error min-imization problem with an objective function that consistsof two parts: physical error, and sparse trajectory error.Using a discretized approximation of the acceleration, theunreduced equation of motion at time step i is given by h − M ( x i − − x i + x i +1 ) = B i +1 p n ,i +1 + f ext , (2)with gravity force f ext . This equation corresponds to our for-ward integration method because the force term evaluationis at the end of the time step. However, we optimize withreduced coordinates z i , where x i = Φ z i , thus, the reducedphysics errors at each time step are C f i ≡ h − Φ T M Φ( z i − − z i + z i +1 ) − Φ T B i +1 p n ,i +1 − Φ T f ext . (3)The desired example trajectory is sparse because it comesfrom an incomplete scan of the surface. Letting vector s i contain the desired point positions at time step i , we canwrite the sparse trajectory error at each time step as C z i ≡ λ ( S Φ z i − s i ) , (4)where the wide sparse selection matrix S extracts the com-ponents of the desired positions by having one non-zeroentry per row. The scalar λ is used to specify the weight ofposition constraints given that the combination of physicsand position constraints are solved in a soft manner.Letting C ( z ) concatenate all physics errors C f on topof all position errors C z , our goal is to find a reducedtrajectory z that minimizes the violation of both, that is,minimizes || C || . We solve this using standard techniques.We do not assemble the Jacobian matrix ∂ C ∂ z directly, but usethe chain rule and keep it in the factored form ∂ C ∂ x ∂ x ∂ z , where ∂ x ∂ z simply contains copies of the basis matrix Φ . The matrix ∂ C f ∂ x has a very simple part that links vertices at different timesteps through the acceleration term, and a more complexpart where the chain rule must be applied to compute theforce gradient. This would normally include a contributionfrom the parametric material correction, but generally wenote better convergence when we omit it, leaving only thegradient of the nominal material. Because we still have theparametric material correction on the right hand side, weonly change the convergence and not the solution. Thus,this second part sprinkles off diagonal terms into the matrixlinking vertices that are adjacent to a common element. Thematrix ∂ C z ∂ x simply contains copies of the selection matrix S . (a) Tracking process (b) Probability of correspondence Fig. 3: Point cloud tracking and shape similarity matching: (a) the physics-based probabilistic tracking method graduallydeforming the mesh to fit the point cloud; (b) the final maximum correspondence for each point in the point cloud to aselection of 10 surface points, which shows tight localized correspondence between the point cloud and the selected points.While the Jacobian matrix is very large, it is also very sparse.We compute the solution using the CUSP [36] library’ssparse least square conjugate gradient implementation onGPU.We can check for convergence to solution z ∗ by mon-itoring our progress in reducing the violation of physicserrors in Equation 3. Once converged, it is exactly theseviolations that provide the necessary control forces to refineour current parametric material correction model. That is,given an optimized reduced trajectory z , the gentle controlforce is computed as f i +1 = h − M Φ( z i − − z i + z i +1 ) − B i +1 p n ,i +1 − f ext . (5)While it may be desirable to solve for control stresses ateach element, as these are what is required for learninga correction, our approach permits an easier solution thatdirectly provides a control force at each vertex. The sparse reduced space-time optimization needs a reason-able starting trajectory. While the forward simulation withthe current parametric material correction could serve thispurpose, we find it valuable to simulate a trajectory that isalso constrained to follow the desired surface motion. Forthe forward simulation, we solve at each step the equation A ∆ v = h f , (6)where f = Bp n + f ext , with p n being the block vector of stresscorrections at the current time step and f ext the externalgravity force, and A = M − hD − h K , where D and K areassembled using the gradient of Equation 1. Many of ourmodels are rigidly attached to the world, and we typicallyremove these degrees of freedom from the system. We canfurther divide the vertices into two groups, (cid:18) A uu A uc A cu A cc (cid:19) (cid:18) ∆ v u ∆ v c (cid:19) = h (cid:18) f u f c (cid:19) , (7)where we use subscripts u for unconstrained and c forconstrained. The second block row can be discarded leavingus a smaller system to solve, namely, A uu ∆ v u = h f u − A uc ∆ v c . (8)Here, ∆ v c at time step i is computed by a second ordercentral finite difference, h − ( x i − − x i + x i +1 ) . We solvethese large sparse systems using PARDISO [37], [38]. For real material observations, the real world data can beincomplete (i.e., only partial surface scans), can be noisy,and can be in the form of unparameterized point clouds.Thus, the exact desired point positions s i at each time step i needed in Equation 4 are no longer available. We mustmodify the trajectory error in the space time optimization toadaptively match nodes of the tetrahedral mesh when theyare close to the surface scans. We do this for every iterationfor finding z ∗ in a manner which is inspired by the physics-based probabilistic tracking method proposed by [39] andextended in [14].For a frame consisting of N points at a given timeinstant (time step i ), we denote point coordinates in thepoint cloud by c n for n = 1 ...N , and vertex positions inthe surface mesh by s k for k = 1 ...K . Let the probabilityof correspondence between the point cloud and the meshvertices be p kn . Assuming that c n is normally distributedaround s k as c n ∼ N ( s k , Σ k ) with an isotropic covariancematrix Σ k = σ I , then we compute the probability thatnodal value s k of the surface mesh corresponds to theobservation c n as p kn = 1 (cid:112) (2 π ) (cid:107) Σ k (cid:107) exp (cid:18) −
12 ( c n − s k ) T Σ − k ( c n − s k ) (cid:19) . Note that parameter σ is set based on the accuracy of thescanner. We only assign point c n to s k if there is a largeprobability (i.e., we truncate p kn to zero if the value isbelow a threshold). Fig. 3 shows how the mesh graduallyconforms to point cloud data, along with final probabilitiesof correspondence between the point cloud and selectedexample vertices.Now, we can reformulate the trajectory constraints ofEquation 4 as a weighed distance between each vertex andits corresponding patch of point clouds as C z i ≡ λ (cid:88) k,n p kn Σ − k ( c n − Φ k z i ) , (9)where Φ k gives the position of surface point k from reducedcoordinates z i , i.e., s k = Φ k z i . The Jacobian matrix requiredfor space-time optimization of this modified position erroris straightforward to compute.During space-time optimization, the correspondenceprobability between point cloud and mesh surface is up-dated, which allows the mesh nodes to move freely across (a) (b) (c) (d) Fig. 4: Von Mises stress visualization for comparison of re-construction results from given nodal elastic forces. The potholder with neo-Hookean material is deformed due to ex-ternal load. (a) Ground truth. (b) RBF based reconstructionresult with 40 kernels. (c) RBF based reconstruction resultwith 80 kernels. (d) Least squares reconstruction result.on the point cloud surfaces. This manner of building param-eterized surface trajectories effectively uses the point clouddata as a soft constraint on the nominal model, and permitsreasonable estimate of surface node trajectories.We use a very similar approach for building the initialstarting point for space time optimization from real data.The tracking procedure in Fig. 3 clearly demonstrates thatthe mesh can well track the point cloud data even in thepresence of large discrepancies.
ARAMETRIC M ATERIAL C ORRECTION
The strain and strain-rate trajectory from the space-time op-timization, combined with the vertex control forces, cannotbe directly used for function fitting. First of all, the controlforces are defined on each vertex, which is highly relatedwith the topology. Moreover, the gentle control forces maynot be entirely self-consistent (i.e., an element might needdifferent forces to correct a given state of strain and strain-rate at different parts of the trajectory). Finally, the amountof data derived from optimization is quite limited. In thissection, we describe how we use two complementary para-metric representations to fit the model.
Section 4.1 introduced the equation B ˆ p = f , which re-lates element stresses to nodal forces. We can compute thestresses as an underdetermined least squares problem usingLSQR [40], which finds ˆ p ∗ that minimizes (cid:107) ˆ p (cid:107) subject to B T B ˆ p = B T f . However, in tests with ground truth data weobserve that the stresses found in this way do not alwaysmatch that well (see Fig. 4). Thus, we use a radial basisfunction representation to regularize the reconstruction ofstress correction through interpolation.From each time step we have a set of principal stretches(and stretch rates) for each element. We assume that wealways have a variety of deformations across elementsand time (i.e., the principal stretches in the data we arefitting are not all the same). With k-means clustering, weselect m stretches to use for interpolation, and define radialbasis functions φ i ( ˆ F ) = || ˆ F − ˆ F i || , for the cluster centers i = 1 , . . . , m . Thus, we can write a smooth strain dependentstress correction as an RBF with m basis functions as ∆ ˆ P = m (cid:88) i =1 w i φ i ( ˆ F ) (10) where w i ∈ R are the weights. For each time step i , wecan assemble a tall matrix R i containing the basis functionsevaluated with ˆ F for each element to write an equation forcomputing interpolated stress reconstructions, ∆ˆ p i = R i w , (11)where the unknown weights are assembled here into a blockvector w . These corrections must explain the gentle forces ofEquation 5, but must also include the current parametricmaterial stress correction ∆ p c , that is, f i = B ( R i w − ∆ p c ) . (12)Without a damping correction, we would solve for theweights with least squares using the above equation at alltime steps. But we include a damping correction througha strain dependent modification of the Rayleigh parameter α , thus, the above equation becomes f i = B ( R i w − ∆ p ) + (cid:88) j ( R ij w α − ∆ α c j ) K j ˙ x , (13)where K j is the contribution of element j to the elementstiffness matrix with its current strain, ∆ α c j is the cur-rent Rayleigh correction for this element, and w α are theweights for computing the new Rayleigh correction with R ij being row j of matrix R i . By least squares we solve for w and w α simultaneously using the data from all time steps. Com-paring with the method of [8], which explicitly interpolatesthe Young’s modulus and Possion ratio using RBF functions,we directly interpolate the relationship between strain andstress. For the damping contribution, however, we use thesame strategy of interpolating physical parameters. The RBF correction provides a smooth parametric materialcorrection, but can also be expensive to compute if manybasis functions are used. Thus, once the iterative materiallearning process has converged, we learn a general neu-ral network representation for the diagonal Piola stresscorrection, N ( ˆ F , ˆ˙ F ) , as seen in Equation 1. The trainingdata are straightforward to compute. We use the learnedRBF parametric model and run many simulations whilerecording principal stretches and stretch rates of elements.During simulation, we also record for each element the totaldiagonal Piola stress correction which is equal to the RBF inEquation 10 plus a dissipation stress contribution computedfrom the ∆ α Rayleigh damping correction.We use a network configuration with two hidden lay-ers as shown in Fig. 5. The network has 6 inputs, theyare principal stretch and its rate, the output is the 3 di-mensional stress correction (total of elastic and dampingstresses). We use 6 nodes in each hidden layer. Havingtested different activation functions with and without batchnormalization, we have settled on using ELU activationfunctions [41] after a batch normalization layer [42] forbetter network performance and more robustness to noise.The evaluation can be seen in Fig. 6, where we used thesame training data for all the network configurations, andthe training data are generated from the first iterationoutput of the turtle example (see supplementary video).
50 100 150 200
Number of epoches T e s t Lo ss ( Log ) -3 SigmoidELUELU+L2BN+ReLUBN+SigmoidBN+ELU
Stretch distance T e s t Lo ss SigmoidELUELU+L2BN+ReLUBN+SigmoidBN+ELU
Fig. 6: Left, the network’s learning rate and test accuracywith different activation function configurations. Right, thenetwork’s test accuracy along deformation scale of test data.We found that the combination of a Batch normalizationlayer + ELU has superior performance and robustness.Fig. 5: Neural network config-uration.As the left plot in Fig. 6shows, most of the acti-vation function configura-tions, except for BN+ReLu,exhibit similar performanceon learning speed and accu-racy when the test data hassimilar deformation scaleas the training data. Herewe use the distance be-tween principal stretch ˆ F and non-deformed principal stretch (1 , , to evaluate thedeformation scale as (cid:107) ˆ F − (1 , , (cid:107) . To show the ability ofa network to extrapolate, we tested the network with muchlarger scales of deformation data, which can be seen in theplot at right of Fig. 6, where the training data resides in theleft side of red dot line, in the range of [0 , . . Beyond thisrange, ELU performs better than the sigmoid function. ESULTS AND D ISCUSSION
In the following sections, we describe experiments that helpreveal what is taking place in each step of the algorithm.To validate the accuracy of our algorithm, we use bothreal captured data and synthetic data generated by forwardsimulations with known material properties.
Space-time optimization is the critical step in the entirepipeline to get training sets from pure kinematic trajectories.For this section we designed three tests using synthetic datato illustrate how our space-time optimization scheme canlead to the convergence of the entire algorithm. To betterreflect real captured data in this evaluation, we also usevirtual scans as input trajectories.For a large scale system or a long trajectory, we mustsolve space-time optimization in reduced space. We testdifferent space-time optimization strategies using the samesynthetic example (turtle with neo-Hookean material) andcompare how well they converge. The results are shown inFig. 7. Tests using node based position constraints or pointcloud patch based position constraints are distinguished bykeyword ‘Node’ and ‘Surface’. From the comparison of the
Material Learning Iteration Number T o t a l N ode - w i s e P o s i t i on E rr o r ( m ) x Traj-Full-NodeTraj-Reduce120-NodeTraj-Reduce60-NodeTraj-Reduce30-NodeTraj-Reduce15-NodePC-Reduce60-Surface
Fig. 7: Convergence comparison for learning a material cor-rection with different space-time optimization strategies. In-put is either synthetic surface trajectories (‘Traj’) or syntheticpoint cloud sequence (‘PC’). Keywords ‘Full’ and ‘Reduce µ ’distinguish between full and reduced space optimization,where µ is the number of basis functions. DenseMiddleSparsePoint Cloud
Material Learning Iteration Number T o t a l N ode - w i s e P o s i t i on E rr o r ( m ) Fig. 8: Convergence comparison under different observationconditions. A neo-Hookean material with edited compres-sion regions is the learning target. From left to right, 100%,13%, and 6% of the surface nodes are assumed visible.Red dots indiacte the exact location of visible nodes. Therightmost bar imitates a real application where the surfaceinformation is represented as point cloud data. The errornorms in nodal positions are plotted on the right withdifferent colors. They all converge sufficiently well.first 5 curves in Fig. 7, node-wise position constraints leadto large control forces during the first several learning iter-ations, when the nominal material is far from ground truth.Consequently, this introduces some overshooting. Throughreduced space-time optimization, these inconveniently largeforce residuals are smoothly distributed throughout theentire object’s domain. We observe good convergence asseen in the last curve in Fig. 7 for learning from point cloudtrajectories.Our algorithm can handle noisy and sparse observa-tions.We tested it on the same synthetic examples withdifferent observation conditions. The learning target is aneo-Hookean material with edited compression regions.Notice that in Fig. 8, the number of observation points inthe first three cases drops sharply. For all these experiments,our algorithm can converge to a sufficiently accurate solu-tion. We also used virtual scans to imitate a real capturesituation, where point cloud data cover the whole surface.Our patch based position constraint adopted in the last case(Section 5.3) performs better than the sparse observationcase and only mildly worse than the accurate full surfaceobservation case.
TABLE 1: Statistics measured for different test cases. From left to right, the test object, the ground-truth constitutive materialmodel, Young’s modulus E G , Poisson’s ratio ν G , Rayleigh damping parameters α G and α G for ground-truth material,nominal material type, the Young’s modulus E N , Poisson ratio ν N , Rayleigh damping parameters α N and α N usedfor nominal model, and maximum position error for testing data reconstruction (for real data this term is the maximumposition error for training data reconstruction). All the Young modulus values are in MPa and all the maximum positionerrors are measured using percentage of object size. All examples use a time step of 0.001 seconds.The letter U denotesuser designed, letter R denotes real materials. Case Material (GT) E G ν G α G α G Material (N) E N ν N α N α N err L Turtle neo-Hookean 2e4 0.43 0.0 0.0 Corotation 3.5e4 0.43 0.0 0.0 2.8Dragon StVK 1e5 0.40 0.0 0.0 Corotation 1.2e5 0.40 0.0 0.0 0.4Sphere1 neo-Hookean(tension) U U 0 0.2 Corotation 1.2e5 0.43 0.0 0.2 4.1Sphere2 neo-Hookean(compression) U U 0 0.2 Corotation 1.2e5 0.43 0.0 0.2 6.0Bar(Damping1) Corotation +strain-dep damping 4e4 0.43 U 0.2 Corotation 4e4 0.43 0.0 0.2 0.1Bar(Damping2) Corotation +strain-dep damping 4e4 0.43 U 0.2 Corotation 4e4 0.43 0.0 0.2 0.8Pot Holder Real Material R R R R Corotation 2.4e6 0.43 0.0 0.0 3.0/7.0Hanger Real Material R R R R Corotation 4.5e5 0.43 0.001 1.0 4.0Silicon Bar Real Material R R R R Corotation 3.0e5 0.45 0 0.0 6.0Bar(Heterogenous) Corotation 1e5/1e7 0.40 0.0 0.0 Corotation 3e6 0.40 0.0 0.0 0.2/0.5 x = 0.01 = 0.1 = 1 = 100 = 50000 = 500000 = 5000000 Material Learning lteration Number T o t a l N o d e - w i s e P o s i t i o n E rr o r ( m ) Fig. 9: Convergence comparison with different position con-straint enforcement weights λ in Eq. (4). A neo-Hookeanmaterial with edited compression regions is the learningtarget. Seven significantly different values are selected, andthe position error norms along learning iteration are plotted.Note that strong position constraints introduce vibration atearly iterations, while soft ones mildly sacrifice accuracy.As described in Equation 4, λ controls the enforcementstrength of position constraints, which influences the con-trol force. We chose different λ values in the wide range[0.01, 500000], and simulated the same example (turtle withneo-Hookean material as learning target, all surface nodesassumed visible). From Fig. 9, we observe that larger λ introduce some vibration or even overshooting at the earlyoptimization stages. The unbanlanced ratio between forceresidual constraint and position constraint causes a largeforce deviation in the internal points, which consequentlydecreases the quality of generated training data. In contrast, the curves corresponding to smaller λ are much smootheralong the optimization iteration axis. From a global perspec-tive, different λ will not necessarily produce significantlydifferent final results, even though there are still subtleaccuracy drops for smaller λ cases as illustrated in zoom-in of Fig. 9. We expect it should be possible to adaptivelytune λ between iterations to speed up convergence. For realcaptured data, λ is chosen to be small, depending on theconfidence in the observations. To validate the generality of our parametric material modelestimation algorithm, we test its ability to learn a variety ofdifferent materials using a co-rotational model for the nom-inal material. Ground truth trajectories are either generatedin the VEGAFem [43] library using classical hyperelastic ma-terial model, or they are user defined elasticity and dampingmodels, or captured by Kinect sensors. Table 1 shows thestatistics of all our testing cases. Qualitative results for thesecases can be seen in the supplementary video. Each case isdiscussed below, while reconstruction errors are listed in thelast column of Table 1.
The turtle (see Fig. 1) is made of neo-Hookean material; thedragon is made of StVK material. We use two deliberatelydesigned test trajectories to validate our learning result. Thefirst test has a similar deformation scale as the trainingtrajectory, while the second test has a much larger rangeof deformation. Table 1 and the supplementary video showthat the learning result can reproduce similar deformation F un c t i on V a l ue f '(x)g'(x)h'(x)f '(x)g'(x)h'(x) F un c t i on V a l ue (a) Learning error for Neo-Hookean material with edited compressed region(b) Learning error for Neo-Hookean material with edited tension region Fig. 10: Visualization of the learning result accuracy for userdesigned elastic material examples [2]. At left the plot showsedited compression and tension regions of a Neo-Hookeanmaterial. At right, the color illustrates principal stress errordistribution of our learned parametric material correction toa corotational model.with high accuracy; the results for different deformationsalso demonstrate low error. Vibration differences can onlybe observed towards the end of a sequence, and as such,can largely be attributed to error accumulation.
Our algorithm can be also extended to model user de-signed nonlinear elasticity and damping material models.Our third example is a soft solid sphere whose top partis keyframed in an up-down motion (Fig. 10). Followingthe method of [2], the internal elastic forces and tangentstiffness matrices are formulated in a polynomial space ofprincipal stretch. Customized materials are designed byediting a single stress-strain curve using a spline interface.We tested two different designed nonlinear materials whichedit tension and compression starting from the a neo-Hookean material as shown in Fig. 10. The supplementaryvideo shows indistinguishable simulation trajectories, whileFig. 10 shows small stress reconstruction errors produced byour parametric material correction of a nominal corotationalmodel. We also compared with the parameter fitting basedmethod of [14] and a modification of the method whichreplaces the simple corotated elasticity model with a neo-Hookean model. The reported side by side comparison inFig. 11 and the supplementary video both demonstrate thatour method is superior to this parameter fitting algorithm,especially when the default model is simple.Although linear viscous damping is widely used inthe computer graphics community, this only constitutes asmall subset of all viscous damping models. Under manycircumstances, the damping matrix C can depend nonlin-early on both deformation and velocity. To validate theaccommodation of our material model estimation algorithmfor damping compensation, Fig. 12 shows our tests on astrain-dependent damping model. We start from a Rayleighdamping model, and substitute the original constant stiff-ness damping coefficient α with a polynomial function of (a) Target[Xu et al. 2015] (b) [Wang et al. 2015]Corotated model (c) [Wang et al. 2015]neo-Hookean model (d) Ours Fig. 11: Comparison with [14] on a user designed materiallearning example. From left to right, (a) target trajectorygenerated using neo-Hookean material with modified com-pressed region using the spline interface from [2]; (b) us-ing [14]; (c) using an enhanced method for [14] that replacesthe corotated model with the neo-Hookean model; and (d)using our algorithm. (a) Deforma�on dependent damping model 1(b) Deforma�on dependent damping model 2
Ground TruthRBF
Ground TruthRBF
Fig. 12: Visualization of the learning result accuracy for userdesigned deformation dependent damping models. Red isused to show ground truth, while blue shows simulatedposes with a learned parametric material correction.the first principal stretch λ . The function of deformationdependent α is controlled using a spline tool.Since our parametric material correction model is inde-pendent of topology, it can be easily transferred to other sim-ulation scenarios. The neo-Hookean material with modifiedtension region model which we learned from ball examplecan be transferred, for example, to a chubby bunny. Ascan be gleaned from the video and Fig. 13(a), the bunnybelly vibrates in a lively fashion during jumps. We alsotransfer the two deformation dependent damping models ofFig. 12 to different leaves of a taro plant to assign separateproperties for young and old leaves as shown in the videoand Fig. 13(b). We first validate our algorithm on a silicon bar example.The object is made by casting an elastomer material oftype Silicon-601 into a bar-shaped mold. The target ma-terial properties depend on the amount of added curingagent,which has no default nominal value. Moreover, thenumerous small air bubbles seen in Fig. 14(d) add morevariance to the material properties. Thus, in the presenceof these bubbles, it is worthwhile to allow an algorithmlike ours to infer the material properties. As illustrated inFig. 14(a) through (c), the object is fixed at one end, andthree different external loads (100 g, 200 g, and 300 g) areadded at the free end. We released the load and recorded itsvibration. The trajectory of the 300 g case is used as trainingdata, while the other two are used for testing. We learn Fig. 13: Top: A chubby jumping bunny using the learned parametric material of the ball example designed material withmodified tension region. Second row: A taro plant responds to user interactions. The left young leaf has damping propertylearned from the top deformation dependent damping model of Fig. 12; while the middle older leaf has damping propertylearned from the second deformation dependent damping model.Fig. 14: Silicon bar examples: (a) g external loading caseis used as training data, (b) g external loading case, and(c) g external loading case are used as test data; (d)irregular bubbles can be observed in a magnified view.the material using both our method and and the methodof [14]. We encourage the readers to watch the side by sidecomparison in the accompanying video.We also validated and compared our algorithm with [14]for the real silicon pot holder and hanger examples, whichare shown in Figs. 15 and 16, respectively. The capturedtrajectories seen in (a) of both figures are used to obtainmaterial corrections. The raw point cloud data are fused,and severe outliers are removed using Artec Studio. Theresults are validated through external loading tests. Morespecifically, we fix the objects on one end, either horizon-tally or vertically, and attach different weights on the otherend. The external weights are suddenly released and thevibrations of the soft object are simulated and comparedwith the ground truth. A side by side comparison canbe seen in the accompanying video. For the hanger andsilicon bar example, we observed that the original constantmass damping coefficient α must be substituted with apolynomial function of principal stretch λ . The algorithm proposed in this paper can also be usedfor material coarsening. In Fig. 17, a high resolution bar (8 × × is composed of two different constitutivematerials, with Young modulus values of 1e5 and 1e7,respectively. The two materials are composited in a layerby layer manner, represented by the light and dark greencolors. The low resolution mesh is the result of coarseningby factor 2 along three axis directions. Two principal defor- (a) (b)20g 500g Fig. 15: Real Material Fitting: (a) Two captured trajectoriesand the corresponding tracking result. (b) Static loadingtests: a silicon pot holder is bent and pulled by externalweights; the holder is fixed at one end horizontally andvertically.Fig. 16: Real Material Fitting: (a) One captured trajectoriesand the corresponding tracking result. (b) Static loadingtests: a silicon hanger is pulled by external weights; thehanger is fixed at top vertically. TABLE 2: Performance statistics measured for different testing cases. Listed from left to right are the test object, numberof vertices, number of tet elements, number of frames for training data, number of reduced modes, number of learningiterations, number of RBF kernels, and total computation time for material learning. The computation times are in hours.
Case
Fig. 17: Material coarsening. The green bar shows the finemesh with a layered material distribution; the purple bar isthe corresponding coarsened mesh with homogeneous ma-terial distribution. Bend (left) and twist (right) deformationtrajectories of fine mesh are used as training data, and thepurple bars are the reconstruction result after learning.mation modes (bend and twist) are used as training data.The equivalent coarsened material property found by ouralgorithm can produce very similar motion to the originalhigh resolution heterogeneous model.
In our algorithm, we use two parametric material repre-sentations. The RBF based representation is used duringiterative material learning while the NN representation istrained after the material correction is achieved. We followstandard practices in training our networks, computingscaling factors for the inputs and outputs based on thetraining data so that both inputs and outputs have zeromean and unit variance. We randomly permute the orderof the samples across time to improve training.Fig. 19 shows the generality of both types of parametricrepresentations. We choose NN, RBF, and RBF augmentedwith a low-order polynomial as candidates. Each row ofFig.19 corresponds to one scenario. The top row has largertest data distribution (blue dots in (a)) than training data(red dots in (a)). All three representations perform well inthe range where the deformations are covered by trainingdata. The performance of both NN and RBF with poly-nomial drops severely when the deformation is out of therange of the training data. In the bottom row, when a sparseset of additional training data is included (pink dots in (e)), (a)(b)
Error (10 -3 ) N u m be r ( l og10 ) Iteration 1Iteration 3Iteration 5Iteration 7Iteration 9Iteration 11
Error (10 -3 ) N u m be r ( l og10 ) Iteration 1Iteration 3Iteration 5Iteration 7Iteration 9Iteration 11
Fig. 18: Comparison of convergence stability showing his-tograms of error sizes with increasing learning iterations, (a)using RBF based representation, the number of position er-rors (forward simulation with material correction comparedto ground truth) in all bins decrease rapidly, except for thesmallest errors (4e-3 and smaller), (b) when using a neuralnetwork representation in the learning method, the errorreduction is slower.the performance for all representations improves greatly.The generality of the different representations is highlydependent on data, especially for complex highly nonlineartarget materials such as the neo-Hookean material withedited tension region that we used here. Since the material isedited in a piece-wise manner, the training data must coverthe entire working space in order to recover the materialfaithfully.Fig. 20 shows that NN training using data that comesfrom simulation sequences leads to unfaithful stress recon-structions as visualized by their isosurfaces in comparisonto ground truth. In contrast, the isosurfaces visualized forthe RBF based representation are close to the required (a) (b) (c) (d)(e) (f) (g) (h) Fig. 19: Generality comparison of parametric representations. Each row corresponds to one specific training dataconfiguration. The first column illustrates the distribution of training (red dots) and test data (blue dots). The pink dots infirst column represent additional sparse training data. The second, third and fourth columns are the result of NN, RBF, andRBF with polynomial’s result respectively, where the Matlab jet colour map is used to indicate the magnitude of the error. (a) (b) (c) (d)
Fig. 20: Comparison on function learning result between two parametric material models. The plots show isosurfaces of theparametric material model for different components (each row shows one component of the principal stress). (a) groundtruth, (b) RBF based representation, (c) NN direct learning result, and (d) NN learning result for fitting the RBF.ground truth correction. Because we can generate muchmore training data through evaluation of our RBF rep-resentation, this leads to a faithful neural network basedrepresentation which is also less expensive to compute.Fig. 18 shows statistics on the distribution of position errors(comparing forward simulation with material correction toground truth) with increasing learning iterations. Using the RBF based representation, the number of larger positionerrors decrease faster with fewer iterations than the errorstabulated for an iterative learning process that uses a net-work model throughout. We measured the computational cost for each critical step ona 10-core 3.0 GHz Intel i7-6950X desktop. The performancefor space-time optimization, listed in Table 2, correlates withthe number of tetrahedral elements, the number of framesin the motion trajectory, the number of RBF kernels, andnumber of selected reduced basis. We also performed aquantitative comparison regarding the computation time forRBF and NN parametric representations. In our test, theRBF representation had kernels, while the NN had thestructure illustrated in Fig. 5. The amount of testing datavaried in the range [200 , with steps of . Accordingto our statistics, the computation time of NN evaluation isconstantly times faster than RBF evaluation under thisconfiguration. ONCLUSION AND F UTURE W ORK
We have presented a new method for estimating nonlinearconstitutive models from trajectories of surface data. Thekey insight is to have a parametric material correction modellearn the error of the elastic and damping properties of anominal material. A framework for gradually learning thiscorrection from only kinematic data is described. We havedemonstrated our method with several examples, illustrat-ing the ability of our approach to learn classical materialmodels, user designed materials (cartoon physics), and realworld captured data.The desire to work with realistic constitutive modelswhen simulating complex motion has been shared for a longwhile by researchers from many fields, not just computergraphics. The possibility of employing machine learningtechnology towards such a goal is tantalizing. Although ourcurrent framework only employs neural network technol-ogy as a compression tool, we solved an interesting andtimely problem which is critical for all machine learningalgorithms for generating annotated data automatically. Webelieve our present work is an important step in that direc-tion.There are a number of interesting ideas to explore in fu-ture research. First, we note that extending our approach toaccommodate a variety of numerical integration techniqueswould help in avoiding or reducing step size dependentnumerical damping effects in our results. Second, we onlyaddress heterogeneous materials in the case of numericalcoarsening to a homogeneous material. There are interestingextensions that can be considered for dealing with varyingmaterial properties across a model, for instance, by addinga latent material parameter to our representation. Finally,there are still a variety of potential damping effects that wecannot capture with our approach. The models we estimatedo not account for any hysteresis in the damping model,while this can be common in real materials, as can alsobe the presence of plastic deformation. Capturing a largervariety of complex plastic and damping behaviors is indeeda very interesting avenue for future work. R EFERENCES [1] H. Xu and J. Barbiˇc, “Example-based damping design,”
ACMTrans. Graph. , vol. 36, no. 4, pp. 53:1–53:14, Jul. 2017. [Online].Available: http://doi.acm.org/10.1145/3072959.3073631 [2] H. Xu, F. Sin, Y. Zhu, and J. Barbiˇc, “Nonlinear materialdesign using principal stretches,”
ACM Trans. Graph. , vol. 34,no. 4, pp. 75:1–75:11, Jul. 2015. [Online]. Available: http://doi.acm.org/10.1145/2766917[3] D. K. Pai, K. v. d. Doel, D. L. James, J. Lang, J. E. Lloyd, J. L.Richmond, and S. H. Yau, “Scanning physical interaction behaviorof 3d objects,” in
Proc. of SIGGRAPH , 2001, pp. 87–96.[4] J. Schoner, J. Lang, and H.-P. Seidel, “Measurement-based inter-active simulation of viscoelastic solids,”
Computer Graphics Forum ,vol. 23, no. 3, pp. 547–556, 2004.[5] M. Becker and M. Teschner, “Robust and efficient estimation ofelasticity parameters using the linear finite element method,” in
Proc. Simulation und Visualization , 2007, pp. 15–28.[6] H. Wang, J. F. O’Brien, and R. Ramamoorthi, “Data-driven elasticmodels for cloth: Modeling and measurement,”
ACM Trans. onGraphics , vol. 30, no. 4, pp. 71:1–71:12, 2011.[7] E. Miguel, D. Bradley, B. Thomaszewski, B. Bickel, W. Matusik,M. A. Otaduy, and S. Marschner, “Data-driven estimation of clothsimulation models,”
Computer Graphics Forum , vol. 31, no. 2, pp.519–528, 2012.[8] B. Bickel, M. B¨acher, M. A. Otaduy, W. Matusik, H. Pfister, andM. Gross, “Capture and modeling of non-linear heterogeneous softtissue,”
ACM Trans. on Graphics , vol. 28, no. 3, pp. 89:1–89:9, 2009.[9] E. Miguel, R. Tamstorf, D. Bradley, S. C. Schvartzman,B. Thomaszewski, B. Bickel, W. Matusik, S. Marschner, and M. A.Otaduy, “Modeling and estimation of internal friction in cloth,”
ACM Trans. Graph. , vol. 32, no. 6, pp. 212:1–212:10, Nov. 2013.[Online]. Available: http://doi.acm.org/10.1145/2508363.2508389[10] E. Miguel, D. Miraut, and M. A. Otaduy, “Modeling andestimation of energy-based hyperelastic objects,”
ComputerGraphics Forum , vol. 35, no. 2, pp. 385–396, 2016. [Online].Available: http://dx.doi.org/10.1111/cgf.12840[11] K. S. Bhat, C. D. Twigg, J. K. Hodgins, P. K. Khosla, Z. Popovi´c,and S. M. Seitz, “Estimating cloth simulation parameters fromvideo,” in
Proc. ACM SIGGRAPH/Eurographics Symp. on ComputerAnimation , 2003, pp. 37–51.[12] S. Yang, J. Liang, and M. C. Lin, “Learning-based cloth materialrecovery from video,” in
The IEEE International Conference onComputer Vision (ICCV) , Oct 2017.[13] A. Davis, K. L. Bouman, J. G. Chen, M. Rubinstein,O. B ¨uy ¨uk¨ozt ¨urk, F. Durand, and W. T. Freeman, “Visual vibrome-try: Estimating material properties from small motions in video,”
IEEE Trans. Pattern Analysis & Machine Intelligence , vol. 39, no. 4,pp. 732–745, 2017.[14] B. Wang, L. Wu, K. Yin, U. Ascher, L. Liu, and H. Huang,“Deformation capture and modeling of soft objects,”
ACM Trans.Graph. , vol. 34, no. 4, pp. 94:1–94:12, Jul. 2015. [Online]. Available:http://doi.acm.org/10.1145/2766911[15] S. Martin, B. Thomaszewski, E. Grinspun, and M. Gross,“Example-based elastic materials,”
ACM Trans. Graph. , vol. 30,no. 4, pp. 72:1–72:8, 2011.[16] C. Schumacher, B. Thomaszewski, S. Coros, S. Martin, R. Sumner,and M. Gross, “Efficient simulation of example-based materials,”in
Proc. ACM SIGGRAPH/Eurographics Symp. on Computer Anima-tion , 2012.[17] D. Chen, D. I. W. Levin, S. Sueda, and W. Matusik, “Data-drivenfinite elements for geometry and material design,”
ACM Trans.Graph. , vol. 34, no. 4, pp. 74:1–74:10, Jul. 2015. [Online]. Available:http://doi.acm.org/10.1145/2766889[18] J. Panetta, Q. Zhou, L. Malomo, N. Pietroni, P. Cignoni, andD. Zorin, “Elastic textures for additive fabrication,”
ACM Trans.Graph. , vol. 34, no. 4, pp. 135:1–135:12, Jul. 2015. [Online].Available: http://doi.acm.org/10.1145/2766937[19] D. Chen, D. I. W. Levin, W. Matusik, and D. M. Kaufman,“Dynamics-aware numerical coarsening for fabrication design,”
ACM Trans. Graph. , vol. 36, no. 4, pp. 84:1–84:15, Jul. 2017.[Online]. Available: http://doi.acm.org/10.1145/3072959.3073669[20] L. Kharevych, P. Mullen, H. Owhadi, and M. Desbrun,“Numerical coarsening of inhomogeneous elastic materials,”
ACM Trans. Graph. , vol. 28, no. 3, pp. 51:1–51:8, Jul. 2009. [Online].Available: http://doi.acm.org/10.1145/1531326.1531357[21] M. Nesme, P. G. Kry, L. Jeˇr´abkov´a, and F. Faure, “Preservingtopology and elasticity for embedded deformable models,”in
ACM SIGGRAPH 2009 Papers , ser. SIGGRAPH ’09. NewYork, NY, USA: ACM, 2009, pp. 52:1–52:9. [Online]. Available:http://doi.acm.org/10.1145/1576246.1531358 [22] R. Torres, A. Rodr´ıguez, J. M. Espadero, and M. A. Otaduy, “High-resolution interaction with corotational coarsening models,” ACMTrans. Graph. , vol. 35, no. 6, pp. 211:1–211:11, Nov. 2016. [Online].Available: http://doi.acm.org/10.1145/2980179.2982414[23] R. M. Snchez-Banderas and M. A. Otaduy, “Strain rate dissipationfor elastic deformations,”
Computer Graphics Forum , vol. 37, no. 8,pp. 161–170, 2018.[24] J. Ghaboussi, J. Garrett Jr, and X. Wu, “Knowledge-based modelingof material behavior with neural networks,”
Journal of EngineeringMechanics , vol. 117, no. 1, pp. 132–153, 1991.[25] J. Ghaboussi, D. A. Pecknold, M. Zhang, and R. M. Haj-Ali, “Au-toprogressive training of neural network constitutive models,”
International Journal for Numerical Methods in Engineering , vol. 42,no. 1, pp. 105–126, 1998.[26] S. Jung and J. Ghaboussi, “Neural network constitutive model forrate-dependent materials,”
Computers & Structures , vol. 84, no. 15,pp. 955 – 963, 2006.[27] D. Stefanos and P. Gyan, “On neural network constitutive modelsfor geomaterials,”
Journal of Civil Engineering Research , vol. 5, no. 5,pp. 106–113, 2015.[28] L. Liang, L. Minliang, and S. Wei, “A deep learning approach to es-timate chemically-treated collagenous tissue nonlinear anisotropicstress-strain responses from microscopy images,”
Acta Biomateri-alia , vol. 63, no. Supplement C, pp. 227 – 235, 2017.[29] R. Luo, T. Shao, H. Wang, W. Xu, X. Chen, K. Zhou, and Y. Yang,“Nnwarp: Neural network-based nonlinear deformation,”
IEEETrans. Visualization & Computer Graphics , 2018.[30] L. Fulton, V. Modi, D. Duvenaud, D. I. W. Levin, and A. Jacob-son, “Latent-space dynamics for reduced deformable simulation,”
Computer Graphics Forum , 2019.[31] E. Sifakis and J. Barbic, “FEM simulation of 3D deformable solids:A practitioner’s guide to theory, discretization and model reduc-tion,” in
ACM SIGGRAPH 2012 Courses , ser. SIGGRAPH ’12, 2012,pp. 20:1–20:50.[32] G. Irving, J. Teran, and R. Fedkiw, “Invertible finite elementsfor robust simulation of large deformation,” in
Proceedings ofthe 2004 ACM SIGGRAPH/Eurographics Symposium on ComputerAnimation , ser. SCA ’04. Aire-la-Ville, Switzerland, Switzerland:Eurographics Association, 2004, pp. 131–140. [Online]. Available:http://dx.doi.org/10.1145/1028523.1028541[33] A. Witkin and M. Kass, “Spacetime constraints,” in
Proceedingsof the 15th Annual Conference on Computer Graphics andInteractive Techniques , ser. SIGGRAPH ’88. New York, NY,USA: ACM, 1988, pp. 159–168. [Online]. Available: http://doi.acm.org/10.1145/54852.378507[34] J. Barbiˇc, M. da Silva, and J. Popovi´c, “Deformable objectanimation using reduced optimal control,”
ACM Trans. Graph. ,vol. 28, no. 3, pp. 53:1–53:9, Jul. 2009. [Online]. Available:http://doi.acm.org/10.1145/1531326.1531359[35] C. Schulz, C. von Tycowicz, H.-P. Seidel, and K. Hildebrandt, “An-imating deformable objects using sparse spacetime constraints,”
ACM Trans. Graph. , vol. 33, no. 4, pp. 109:1–109:10, Jul. 2014.[Online]. Available: http://doi.acm.org/10.1145/2601097.2601156[36] S. Dalton, N. Bell, L. Olson, and M. Garland, “Cusp:Generic parallel algorithms for sparse matrix and graphcomputations,” 2014, version 0.5.0. [Online]. Available: http://cusplibrary.github.io/[37] C. G. Petra, O. Schenk, M. Lubin, and K. G¨artner, “An augmentedincomplete factorization approach for computing the Schur com-plement in stochastic optimization,”
SIAM Journal on ScientificComputing , vol. 36, no. 2, pp. C139–C162, 2014.[38] C. G. Petra, O. Schenk, and M. Anitescu, “Real-time stochasticoptimization of complex energy systems on high-performancecomputers,”
IEEE Computing in Science & Engineering , vol. 16, no. 5,pp. 32–42, 2014.[39] J. Schulman, A. Lee, J. Ho, and P. Abbeel, “Tracking deformableobjects with point clouds,” in
Proc. IEEE Int. Conf. on Robotics &Automation , 2013.[40] C. C. Paige and M. A. Saunders, “LSQR: An algorithm for sparselinear equations and sparse least squares,”
ACM Trans. Math.Softw. , vol. 8, no. 1, pp. 43–71, Mar. 1982. [Online]. Available:http://doi.acm.org/10.1145/355984.355989[41] D.-A. Clevert, T. Unterthiner, and S. Hochreiter, “Fast and accu-rate deep network learning by exponential linear units (elus),”
Computer Science , 2015.[42] S. Ioffe and C. Szegedy, “Batch normalization: Accelerating deepnetwork training by reducing internal covariate shift,” in