Adv Search
Home | Accepted | Article In Press | Current Issue | Archive | Special Issues | Collections | Featured Articles | Statistics

2021,  3 (2):   156 - 170   Published Date:2021-4-20

DOI: 10.1016/j.vrih.2021.01.002
1 Introduction2 Related work2.1 Fast deformable simulation2.2 Inhomogeneous material simulation3 Rationale and overview4 Nonlinear homogenization4.1 Nonlinear deformation modes4.2 Anisotropic nonlinear material model4.3 Optimization-based homogenization4.4 Implementation details5 Results6 Conclusions


Fast simulation techniques are strongly favored in computer graphics, especially for the nonlinear inhomogeneous elastic materials. The homogenization theory is a perfect match to simulate inhomogeneous deformable objects with its coarse discretization, as it reveals how to extract information at a fine scale and to perform efficient computation with much less DOF. The existing homogenization method is not applicable for ubiquitous nonlinear materials with the limited input deformation displacements.
In this paper, we have proposed a homogenization method for the efficient simulation of nonlinear inhomogeneous elastic materials. Our approach allows for a faithful approximation of fine, heterogeneous nonlinear materials with very coarse discretization. Modal analysis provides the basis of a linear deformation space and modal derivatives extend the space to a nonlinear regime; based on this, we exploited modal derivatives as the input characteristic deformations for homogenization. We also present a simple elastic material model that is nonlinear and anisotropic to represent the homogenized materials. The nonlinearity of material deformations can be represented properly with this model. The material properties for the coarsened model were solved via a constrained optimization that minimizes the weighted sum of the strain energy deviations for all input deformation modes. An arbitrary number of bases can be used as inputs for homogenization, and greater weights are placed on the more important low-frequency modes.
Based on the experimental results, this study illustrates that the homogenized material properties obtained from our method approximate the original nonlinear material behavior much better than the existing homogenization method with linear displacements, and saves orders of magnitude of computational time.
The proposed homogenization method for nonlinear inhomogeneous elastic materials is capable of capturing the nonlinear dynamics of the original dynamical system well.


1 Introduction
Since Terzopoulos et al.'s seminal paper[1] on elastically deformable models, physical-based deformable simulation has become popular in computer graphics, with tremendous developments being made. Despite the progress that has occurred, most of the work has focused on objects made of a single, homogeneous material. This is because simulation methods with this simplification are more applicable to interactive applications. However, many real-world objects are composed of heterogeneous materials. Simulation of such complex objects using the currently available techniques usually requires high-resolution spatial discretization to resolve the fine-scale heterogeneity. This requirement leads to overwhelming computational costs, thereby making interactive simulation impractical.
Our goal is to efficiently simulate inhomogeneous deformable objects with very coarse discretization, while effectively capturing the physical behavior. Homogenization theory[2] is a perfect match for this purpose, as it reveals exactly how to extract information at a fine scale to perform efficient computation. However, the application of homogenization theory in computer graphics is surprisingly rare. To the best of our knowledge, Kharevych et al. firstly introduce the homogenization method for graphics applications[3]. From the exact matching of potential energies for a set of characteristic displacements, their method achieves a coarse approximation of deformable objects that are composed of inhomogeneous linear elastic materials. The coarsened material properties allow for real-time simulation that captures the correct dynamic behavior of the original materials. Their method is limited to linear elasticity; therefore, it is not applicable for ubiquitous nonlinear materials. In this paper, we address this problem and propose a novel homogenization method for nonlinear inhomogeneous elastic materials.
We employed the same strategy as in the method of Kharevych et al., obtaining the coarsened material properties by matching the potential energies for a set of characteristic displacements. The homogenization of nonlinear materials imposes several key challenges compared to linear homogenization. First, the space of nonlinear deformations is significantly larger than the space of infinitesimal linear deformations. Therefore, a set of displacements that are characteristic of typical nonlinear deformations has to be found. Second, a proper anisotropic nonlinear material model needs to be defined for the coarse representation so that the homogenized material properties exhibit sufficient anisotropy and nonlinearity. Finally, the exact matching of potential energies for the chosen displacements requires that the number of equations is equal to the number of unknowns in the material properties to be solved. This requirement is overly restrictive for nonlinear homogenization because nonlinear materials possess a larger family of characteristic deformations, and it is good practice to account for as many of them as possible during homogenization.
To address each of these problems, the respective contributions are presented in this paper. We exploited the idea of using deformation modes from modal analysis as the typical deformations for homogenization. Although the modal basis has been explored extensively for reduced simulation in previous methods[4], using the inherent properties of the modal basis for homogenization is new. Nonlinear deformations are beyond the span of the linear modal basis. We constructed the deformation basis employing the modal derivatives technique proposed by Barbič and James[5], and thus, nonlinear deformations are covered by the basis.
Our second contribution to the field is a simple anisotropic nonlinear material model. Based on the isotropic St. Venant-Kirchhoff model, which is commonly used in graphics, we defined its anisotropic extension. This anisotropic material model was used in our coarse representation and expressed sufficient nonlinearity for graphics simulations.
Furthermore, unlike the previous method that solves for the coarsened material properties with exact matching of the potential energies, we formulated the homogenization process as an optimization. The weighted sum of the strain energy deviations for the input deformation modes was minimized. We also considered the frequencies of the deformation modes and put greater weight on the more important low-frequency modes. With this optimization strategy, an arbitrary number of bases can be used as input data, and therefore, the restriction on the number of displacements can be overcome.
2 Related work
A complete review of the vast number of methods for the simulation of physical deformation is beyond the scope of this paper, and we refer the reader to the excellent survey by Nealen et al.[6]. In this section, we focus on previous research closest to ours on fast simulation techniques and the simulation of inhomogeneous materials.
2.1 Fast deformable simulation
Fast simulation techniques are strongly favored in computer graphics because real-time and interactive frame rates are crucial to applications such as video games and virtual surgery. Many strategies have been devised for fast simulation. Running simulations with larger time intervals[7] is one strategy that seeks to reduce the overall computation cost by using fewer simulation steps. Other methods[8-10] exploit the simplicity of linear material models while circumventing the artifacts with large deformations by factoring out rotation from the deformation before applying the material model and then rotating it back. These methods are called stiffness warping methods or linear corotated methods.
Multi-resolution approaches are intuitive solutions that adaptively refine the simulation for more intensive deformations and coarsen it otherwise. The homogenization method is one such method; it utilizes a low-resolution mesh to animate the approximate motion of the high-resolution mesh. The advantage of the homogenization method is that it improves the simulation efficiency by using a coarse mesh instead of a fine mesh, which easily reduces the degrees of freedom (DOF) and dimensions entailed in solving linear equations. Because objects typically interact through their surfaces, another way of reducing the DOF in simulations is to express the physical equations at the surface points only[11-12]. However, such methods are generally limited to small deformations. In graphics, domain embedding[13-14] is a widely used technique that achieves fast computation and detailed rendering by embedding high-resolution visual meshes into a coarse deformable simulation.
Modal reduction methods limit the possible deformations in a low-dimensional subspace to achieve efficiency. The reduced space is constructed using a set of deformation bases, typically obtained through eigenvalue analysis of the stiffness matrix or of empirical simulation data. Among the two strategies, methods using eigen-analysis of the stiffness matrix, or so-called modal analysis, are more prevalent. Linear modal analysis provides the deformation basis only for small deformations from the configuration at which the stiffness matrix is evaluated, and is therefore not suitable for the simulation of large deformations. Several solutions have been proposed to resolve this problem. Choi and Ko exploited the idea of stiffness warping for reduced simulation and developed a procedure to express and update the rotation component of deformation in a subspace[15]. Barbič and James generalized linear modal analysis and proposed an enhancement of the linear deformation modes with their directional derivatives[5]. These so-called modal derivatives contain substantial nonlinear content that is sufficient to describe large deformations in practice. The idea of modal reduction is extensively used by subsequent methods to reduce computation complexity[16-19].
Instead of using a modal basis for reduced simulation, in this study, we exploited the inherent properties of modal bases as characteristic deformations and used them for homogenization. To our knowledge, modal analysis has not been used for this purpose.
2.2 Inhomogeneous material simulation
Simulation of deformable objects that are composed of inhomogeneous materials is less explored in computer graphics owing to the complexity of the material structure. Nesme et al. proposed approximating non-uniform stiffness on cubic grids with the spatial average of the elasticity tensor[20]. However, such a simple average does not accurately coarsen the original materials. In another work, Nesme et al. employed a domain embedding technique and computed shape functions for coarse elements, considering the varying materials in the elements based on high-resolution mechanical analysis in the pre-computation step[21]. They also considered the void space inside the elements and preserved the branch structure with element duplication. Unfortunately, their method only applies to linear elasticity. Bickel et al. presented a data-driven approach to simulate nonlinear heterogeneous materials[22]. With a set of stress-strain relationships captured from example deformations, they modeled the material by nonlinear interpolation of the relationships at run-time. As the captured data are restricted to simple deformations with a single contact probe, their method cannot accurately capture the material behavior in complex scenarios. Faure and colleagues introduced material-aware shape functions that efficiently resolve non-uniform stiffness for sparse meshless skinning[23]. These shape functions failed to ideally resolve complex three-dimensional deformations because they only used one scalar stiffness value for computing the shape functions. Kharevych et al. creatively adopted homogenization theory in computer graphics and properly approximated heterogeneous materials on very coarse discretizations with orders of magnitude speedup[3]. Their method is based on linear elasticity and is therefore restricted to infinitesimal deformations of linear materials. In this paper, we have addressed these limitations using a nonlinear homogenization method. Chen et al. proposed building a material model database to avoid the repeated calculation of material parameters in the simulation process[24]. The proposed method only queries the coarse resolution material information in the database corresponding to the high-resolution mesh. This method can simulate more complex nonlinear materials, but it ignores the relationship between the elements. At the same time, the material model depends on the training set and parameter selection process. In their subsequent research, they proposed a material solution scheme that considered the inconsistency of motion frequency and boundary processing problems in the simulation process of different resolutions, which achieved good simulation results[25]. Chen et al. proposed the use of the matrix shape function to replace the traditional scalar shape function to simulate the nonlinear constitutive model of heterogeneous materials[26]. In the pre-computation process, the balance of the geometric continuity and local material rigidity information is considered. In contrast with Chen et al.'s previous method[25], the advantages of this method entail the deformation characteristics of heterogeneous materials without adjusting parameters. In the following year, they used a wavelet transform to propose a shape function that supported local features and improved the boundary problems[27]. However, they use a linear corotated model, and there is still a certain expansion space for the simulation of the nonlinear constitutive model. In this study, we aimed to build a nonlinear and anisotropic material model by exploiting modal derivatives as the nonlinear input displacement to generate a more diverse deformation space.
3 Rationale and overview
In this section, we present a brief overview of our homogenization algorithm. For consistency and clarity, we use SOLID characters to refer to quantities on the fine discretizations and
characters to refer to their coarse counterparts.
We employed FEM discretization for simulations; therefore, our method started from a fine mesh
with varying material properties and obtained an approximation of the material properties on a coarse mesh
whose element count is much smaller. Following Kharevych et al.'s approach[3], we computed the coarse material parameters by matching the potential energies between two scales for a set of characteristic deformations. Instead of enforcing exact matching, we tolerated some deviation in the energies. We defer the explanation of the benefits to the section that follows, and merely describe the methodology here. For any element
on the coarse mesh, several elements,
T q
, on the fine mesh occupy the same space. Hence, the potential energy of the coarse element
should approximately match the sum of the potential energies over the fine elements
T q
for all deformations.
W ( 𝕋 P ) W ( T q )
Computing the coarse potential energy
W ( 𝕋 P )
requires a proper constitutive model defined on the coarse mesh. As inhomogeneous materials generally lead to anisotropic behavior, the material model must be anisotropic. To represent homogenized nonlinear materials, the material model also needs to be nonlinear. We defined a simple material model in this study to achieve these aims.
Kharevych et al. used material deformations subjected to linear tractions on the boundary as the input displacements. In contrast, we employed a nonlinear modal basis to enforce the energy matching. As typical nonlinear deformations were taken into account by our homogenization process, we wished to obtain graceful approximations for general nonlinear deformations. Given a displacement field
on the fine mesh, we could easily downsample it via interpolation and obtain the displacement field
on the coarse mesh. For boundary nodes of the coarse mesh that lie outside of the fine mesh domain, we found the closest fine element to it and used extrapolation. The potential energies for meshes with different resolutions were evaluated using
, respectively.
To find the elements from a fine mesh for each corresponding coarse mesh, we needed to detect the intersection tetrahedral elements. To build up the energy equivalence of these two meshes with different resolutions, we had to know the intersection volume of each fine element and the coarse tetrahedral mesh. We simply chose the Monte Carlo approximation to solve the intersection volume of each tetrahedral element. We sampled M points in the coarse element and determined the corresponding number of fine elements. Then, the intersection volume of each fine element in the coarse element is
  V f i n e = N M V c o a r s e
, where V denotes the volume of the tetrahedral element. Although the simplification incurs some time overhead, our material solving process is a pre-computation process, and the efficiency is acceptable.
4 Nonlinear homogenization
The strength of our method is achieved by combining three contributions: the use of nonlinear deformation modes, an anisotropic nonlinear material model, and an optimization-based homogenization strategy. We now present detailed descriptions of our work.
4.1 Nonlinear deformation modes
Kharevych et al.'s approach and ours share the same basic idea of enforcing potential energy matching between two scales for a set of characteristic deformations, expecting that the coarsened material matches the original heterogeneous one for all possible deformations. It is important to find a set of deformations that are typical enough to represent all deformations, that is, they form the basis of a deformation space.
Deformation basis generation is a difficult open problem, especially for deformations under general forcing. Kharevych et al. used a set of so-called harmonic displacements as typical deformations for homogenization. These displacements were computed by solving a set of static boundary value problems with linear surface tractions prescribed as Neumann boundary conditions.
Inspired by recent methods that used modal analysis to construct a reduced deformation space, we exploited the fact that modal basis represents typical deformations from the frequency perspective. The standard linear modal basis can only express small deformations, and we employed the solution provided by Barbič and James[5], which constructed a set of nonlinear deformation modes by including directional derivatives of linear modals. Their paper provides a detailed description of the procedures to compute modal derivatives and to generate the low-dimensional deformation basis with mass-PCA. Our contribution lies in a new application of modal analysis for homogenization. We took the homogeneous bar model as an example to demonstrate the difference of the displacement (Figure 1). Harmonic displacements merely consist of linear stretch and shear (top); therefore, they encode no information on nonlinear deformations, whereas our nonlinear deformation modes (bottom) contain substantial nonlinearity.
As the homogeneous bar model shows in Figure 1, the top row demonstrates the harmonic displacements of Kharevych et al.'s approach, which merely consists of linear stretch and shear, and the bottom row adopts nonlinear deformation modes by using modal basis, which contains substantial nonlinearity.
4.2 Anisotropic nonlinear material model
As mentioned earlier, the material model for the coarse mesh needs to be anisotropic and nonlinear in order to express complex behaviors due to heterogeneity and nonlinearity. However, the wide variety of material models in the literature, such as the popular Neo-Hookean and Mooney-Rivlin models, are mostly isotropic. Here, we have defined an anisotropic and nonlinear material model based on the simple yet widely used St. Venant-Kirchhoff model.
The St. Venant-Kirchhoff material, or StVK in short, is perhaps the simplest type of material model. It extends the isotropic linear elasticity model
σ = λ ( t r ϵ ) I + 2 μ ϵ
with geometric nonlinearity, replacing the Cauchy strain
with the Green strain
. Its material behavior can be expressed by the relationship between
and the second Piola-Kirchhoff stress
S = λ ( t r E ) I + 2 μ E
. The stress-strain relationship is linear, but the nonlinear strain tensor injects sufficient nonlinearity between displacements and stress for applications in computer graphics.
Our anisotropic material model exploited the simplicity of the StVK model with regard to geometric nonlinearity, and added anisotropy by relating the stress and strain with a rank-4 elasticity tensor
instead of two scalars:
S = C   :   E
denotes the double contraction of tensors. This material model can also be regarded as an extension of the linear model used by Kharevych et al., but with geometric nonlinearity. The density of the potential energy is computed as
Ψ = S   :   E
. As a result, our homogenization algorithm also solved for the elasticity tensor
𝕋 p
on the elements
𝕋 p
of the coarse mesh. The inclusion of geometric nonlinearity in the homogenization process worked well, which we have demonstrated with a comparison of simulations between two scales under large deformations.
4.3 Optimization-based homogenization
Enforcing the exact equivalence of strain energies between each coarse element and its corresponding fine elements was unnecessary because the error introduced by discretization was inevitable. The best we could expect was to minimize the deviation between the two scales. In addition, exact equivalence requires that the number of equations equals the number of unknowns in the material properties to be solved. In the case of solving for a symmetric elasticity tensor
, the number of input displacements must be six. Six displacements might be sufficient to describe typical linear deformations, but it is not for large deformations because of the significantly larger deformation space.
Consequently, we formulated homogenization as an optimization problem with respect to the elasticity tensor
. We further explored the importance of different deformation modes and put greater weight on low-frequency modes. The objective to be minimized for each coarse element
𝕋 p
was the weighted sum of the potential energy deviations over all input deformations. The weight for each deformation mode was determined using the same strategy as Barbič and James[5] to select the basis during mass-PCA. Considering the physical property of the elasticity tensor,
is represented as an invertible symmetric
6 × 6
matrix with non-negative entries. As a result, the optimization with constraints for the coarse element
𝕋 p
m i n 𝕋 p i = 1 r 1 2 ω i ( W ( 𝕌 i ) - T q D T q 𝕋 p 0 | T q 𝕋 p | | 𝕋 p | W ( u i ) ) 2   s . t .            𝕋 p ( k ,   l ) 0 , 1 k ,   l 6            𝕋 p = 𝕋 p T            d e t ( 𝕋 p ) 0
is the number of input deformation modes, and
u i
𝕌 i
are the i-th deformation mode on the fine mesh and the coarse mesh, respectively. The
operator computes the volume of the underlying region. The i-th deformation mode is either a linear mode or a modal derivative, and the weight
ω i
is determined as follows:
ω i = λ 1 λ j , i f u i = Ψ j λ 1 2 λ j λ k , i f u i = Φ j k
is the eigenvalue corresponding to the linear deformation modes, and
stand for the linear modes and modal derivatives.
We used the solver provided by NLopt[28] to solve this optimization, and different coarse elements could be solved in parallel to save computation time.
4.4 Implementation details
Multiple-resolution mesh generation. In our homogenization algorithm, we generated a low-resolution simulation model from the high-resolution model. We prepared meshes with two resolutions as input data. For structured models, we generated boundary aligned meshes manually or utilized the software MeshLab to simplify the fine surface mesh and generate a coarse mesh by remeshing to maintain its quality for the numerical computing in the simulation process. Then, we generated the tetrahedral meshes of the two meshes. For complex models, it is difficult to maintain key structures of the model using the mesh simplification operation and remeshing method.
Figure 2 shows an example of our homogenization procedure setup using a plant as a model. The stiffness distribution of the plant model is visualized using a color map (left), and our homogenization algorithm obtained material properties on a very coarse mesh from material properties on a highly detailed mesh. A plant is divided into three basic parts: the root is the hardest with Young's Modulus 5.0e6, the leaves are the softest part with Young's Modulus 1.0e6, and the Young's Modulus of the stem is 5.0e5.
For the model with complex structures, we utilized the voxelization method[5] to obtain the volumetric mesh of the simulation meshes with two resolutions. However, with this method, the two meshes have a boundary-conforming problem. To solve the inconsistency of the boundaries, we proposed a simple solution. We generated the surface mesh of the voxelization low-resolution mesh, and computed the signed distance to the high-resolution mesh. Then, the vertex of the low-resolution mesh was moved to the high-resolution mesh in the opposite direction to the gradient of the distance field. After the adjustment process, the meshes with these two resolutions achieved higher boundary consistency, which could be satisfied with our algorithm. Finally, we generated meshes with two resolutions, as depicted in Figure 2 (high resolution: middle, low resolution: right). In addition to our simple method, Stuart et al. proposed another higher tightly volumetric mesh from a high-resolution surface mesh[29].
FEM simulation. To validate our algorithm, we used the Vega FEM library[30] to run simulations at different resolutions. As the high-resolution mesh was sufficiently dense, we assumed that the material model for each fine mesh element was isotropic. The global behavior was still anisotropic because the parameters of the elements were different. We employed the StVK model for each fine element in the high-resolution simulations of nonlinear materials, while our anisotropic material model was used in the corresponding coarse simulations. For comparison with Kharevych et al.'s method, we ran simulations of linear materials using corotated methods with isotropic materials at high resolutions and anisotropic materials at coarse resolutions. The backward Euler integration method with a uniform time step of
Δ t = 0.001 s
was used for all simulations.
Rendering. We embedded surface meshes with high-quality details in the simulation meshes for rendering purposes. The surface meshes were deformed during the simulation via interpolation of the simulation mesh positions. The deformed surface meshes were rendered offline to generate the figures presented in this paper.
5 Results
We have demonstrated the power of our method with several carefully designed examples. High-resolution material models are nonlinear and the objects typically undergo large deformations. By comparing with the results of Kharevych et al.'s method, we have shown that the homogenized material properties obtained from our method approximate the original nonlinear material behavior much better than Kharevych et al.'s linear approach. The Poisson ratio is 0.3 for all examples.
In Figure 3, we fixed the top end of an inhomogeneous elastic bar with layered materials and lifted its bottom end horizontally. The bar swings with large stretches and bending after the lifted end is released. A side-by-side comparison of high-resolution and low-resolution simulations indicates that our homogenized material is more faithful to the original heterogeneous material. An inhomogeneous layered bar fixed at the top swings under gravity. The regions in blue (Young's Modulus:2.0e6) are stiffer than those in yellow (Young's Modulus:5.0e5).
Figure 4 is another example of inhomogeneous deformable objects undergoing nonlinear deformations. The elastic bar made of nine layers deforms with coupled twisting and stretching. The Young's Modulus of the stiffer part is 5.0e5 (blue) and that of the softer part is 1.0e5 (yellow). The homogenized model of Kharevych et al.'s method fails to resolve such nonlinearity, and the coarse simulation exhibits a significant difference from the fine simulation. In contrast, coarse simulation using our results is consistent with the high-resolution simulation.
Figure 5 shows a T-shaped structure with Young's Modulus 3.5e6, where the vertical direction is fixed, and the fiber structure with higher stiffness (Young's Modulus: 5.0e7) is embedded in different directions in the horizontal direction. The complicated T-shaped model deforms under gravity. Because of the different embedding fibers, the two horizontal branches show different material properties. We compared our method with that of Kharevych et al., and the different simulation frames are shown in Figure 5.
Figures 6 to 8 demonstrate the comparison between our method and Kharevych et al.'s method based on the quantitative analysis of bar swing, bar twist, and T-shape examples. A vertex with large deformation was selected on the surface mesh of the object, and the displacement curve of the point in the vertical direction was drawn with different resolutions and the material properties obtained by different methods. The deviation between the variation curves corresponding to the high- and low-resolution simulations reflects the error between the results of the homogenization algorithm and the reference value. As can be seen from the figure, the deviation between the curves corresponding to our method is much smaller than that corresponding to Kharevych's method.
Table 1 lists the timing information for the presented examples. The computations were performed on a PC with a quad-core Intel Core i5, 2.8 GHz CPU. Homogenization was conducted in parallel using four CPU cores, while the simulations ran sequentially on a single core. As can be seen from the table, with homogenization as a pre-computation process, we achieved orders of magnitude speedup in the simulation. All the coarse simulations presented in this paper ran at real-time frame rates.
Summary of the results
Example Ef Ec th(min) tf(ms) tc(ms) Speedup
bar swing 13056 702 19.0 221.9 9.2 24.1x
bar twist 13056 702 21.0 215.1 10.5 20.5x
T 44413 584 49.9 846.3 8.1 104.5x
The columns indicate the number of elements on the fine mesh and coarse mesh, the time for homogenization (in minutes), the simulation time per time step (in milliseconds) for high-resolution simulation, coarse simulation, and the speedup ratio of the coarse simulation to high-resolution simulation, respectively.
6 Conclusions
We have presented a homogenization method that can obtain the effective material properties of nonlinear inhomogeneous materials on very coarse discretizations. Coarse simulations with our homogenized materials captured the original material behavior with acceptable accuracy and saved orders of magnitudes of computation time.
In theory, the coarse mesh and the fine mesh yield a relationship between geometry representation and deformation energy; thus, the proposed method is valid for both microstructural and structural models. However, during the experimental process, we found that the algorithm employed in this study still has some limitations. For the model with a complex shape and structure, when the deformation was complex, the deformation effect of our homogenization algorithm on the low-resolution mesh could reflect the stiffness distribution of the material. Although the trend of motion was consistent with that obtained using a high resolution, there were still obvious differences with the high-resolution case in terms of frame comparison. As depicted in Figure 9, a plant swings with an initial deformation. The figure shows obvious differences between the simulation effect of the low-resolution mesh obtained with our homogenization algorithm (row 2) and that of the high-resolution mesh (row 1). The stiffness of the plant simulation object is shown in Figure 2 (left-1). After analysis, it was found that the displacement deviation may be due to the inherent error caused by finite element simulation and the discrete solutions of meshes with different resolutions. Therefore, to verify this possibility, we used the same homogeneous material in the high- and low-resolution meshes of the plant model, as shown in Figure 10. The top row shows the high-resolution simulation, while the bottom row shows the low-resolution simulation. Owing to the mismatching of geometric boundaries and the errors introduced by discretization, even if the two resolution meshes have the same parameters, the simulation results still have differences. In order to eliminate the differences in deformation frequency caused by different discrete meshes, the authors of [25] proposed a practical method.
Our goal was to enhance the simulation efficiency by using a coarse mesh with inhomogeneous materials instead of the original mesh with a high resolution. Our experimental results demonstrated that our method has attained this goal while maintaining a deformable behavior similar to that of the high-resolution model to some extent. Therefore, we created a coarse mesh with a lower resolution for each experiment instead of testing the resolution limitation of the coarse mesh to generate simulation effects similar to those of the original high-resolution mesh.
In our method, we considered the corresponding factors such as the displacement, volume, and meshes in the precomputing process. For the elastic tensor calculation process, we mainly considered the energy equivalence relationship between each coarse element and its corresponding fine elements. It is possible to consider the local stress fields corresponding to the relationship between the two resolution meshes or the displacements of the corresponding vertices as a constrained condition in our future work.
It would also be interesting to explore the fast simulation of inhomogeneous materials undergoing topology changes using homogenization techniques. The material properties of the corresponding elements need to be recomputed at run-time in case of topology changes. Our optimization-based homogenization strategy is efficient, whereas much work remains to be done in order to re-homogenize local elements on the fly without jeopardizing the simulation.



Terzopoulos D, Platt J, Barr A, Fleischer K. Elastically deformable models. ACM SIGGRAPH Computer Graphics, 1987, 21(4): 205–214 DOI:10.1145/37402.37427


Gloria A. Numerical homogenization: survey, new results, and perspectives. ESAIM: Proceedings, 2012, 37: 50–116 DOI:10.1051/proc/201237002


Kharevych L, Mullen P, Owhadi H, Desbrun M. Numerical coarsening of inhomogeneous elastic materials. In: ACM SIGGRAPH 2009 papers on-SIGGRAPH '09. New Orleans, Louisiana, New York, ACM Press, 2009, 28(3): 51:1–51:8 DOI:10.1145/1576246.1531357


Sifakis E, Barbic J. FEM simulation of 3D deformable solids: a practitioner's guide to theory, discretization and model reduction. In: ACM SIGGRAPH 2012 Posters on-SIGGRAPH '12. Los Angeles, California, New York, ACM Press, 2012, 1–50 DOI:10.1145/2343483.2343501


Barbič J, James D L. Real-Time subspace integration for St. Venant-Kirchhoff deformable models. SIGGRAPH '05: ACM SIGGRAPH 2005 Papers, 2005, 982–990 DOI:10.1145/1186822.1073300


Nealen A, Müller M, Keiser R, Boxerman E, Carlson M. Physically based deformable models in computer graphics. Computer Graphics Forum, 2006, 25(4): 809–836 DOI:10.1111/j.1467-8659.2006.01000.x


Gast T F, Schroeder C, Stomakhin A, Jiang C, Teran J M. Optimization integrator for large time steps. IEEE Transactions on Visualization and Computer Graphics, 2015, 21(10): 1103–1115 DOI:10.1109/tvcg.2015.2459687


Müller M, Dorsey J, McMillan L, Jagnow R, Cutler B. Stable real-time deformations. In: Proceedings of the 2002 ACM SIGGRAPH/Eurographics symposium on Computer animation-SCA'02. San Antonio, Texas, New York, ACM Press, 2002, 49–54 DOI:10.1145/545261.545269


Müller M, Gross M. Interactive virtual materials. Proceedings of Graphics Interface. London, 2004, 239–246


Georgii J, Westermann R. Corotated finite elements made fast and stable. Proceedings of the 5th Workshop on Virtual Reality Interactions and Pysical Simulations. VRIPHYS 2008, Grenoble, France, 2008,11–19


James D L, Pai D K. Real time simulation of multizone elastokinematic models. In: Proceedings 2002 IEEE International Conference on Robotics and Automation. Washington, DC, USA, IEEE, 2002, 927–932 DOI:10.1109/robot.2002.1013475


James D L, Pai D K. Multiresolution Green's function methods for interactive simulation of large-scale elastostatic objects. ACM Transactions on Graphics, 2003, 22(1): 47–82 DOI:10.1145/588272.588278


Capell S, Green S, Curless B, Duchamp T, Popović Z. Interactive skeleton-driven dynamic deformations. ACM Transactions on Graphics, 2002, 21(3): 586–593 DOI:10.1145/566654.566622


Teschner M, Heidelberger B, Muller M, Gross M. A versatile and robust model for geometrically complex deformable solids. In: Proceedings Computer Graphics International. Crete, Greece, IEEE, 2004, 312–319 DOI:10.1109/cgi.2004.1309227


Choi M G, Ko H S. Modal warping: real-time simulation of large rotational deformation and manipulation. IEEE Transactions on Visualization and Computer Graphics, 2005, 11(1): 91–101 DOI:10.1109/tvcg.2005.13


Barbič J, Zhao Y. Real-time large-deformation substructuring. Acm Transactions on Graphics, 2011, 30(4):91:1–91:8 DOI:10.1145/1964921.1964986


Hildebrandt K, Schulz C, von Tycowicz C, Polthier K. Interactive spacetime control of deformable objects. ACM Transactions on Graphics, 2012, 31(4): 71 DOI:10.1145/2185520.2185567


Harmon D, Zorin D. Subspace integration with local deformations. ACM Transactions on Graphics, 2013, 32(4): 1–10 DOI:10.1145/2461912.2461922


von Tycowicz C, Schulz C, Seidel H P, Hildebrandt K. An efficient construction of reduced deformable objects. ACM Transactions on Graphics, 2013, 32(6): 213 DOI:10.1145/2508363.2508392


Nesme M, Payan Y, Faure F. Animating shapes at arbitrary resolution with non-uniform stiffness. Proceedings of the 3rd Workshop in Virtual Reality Interaction and Physical Simulation. Aire-la-Ville: Eurographics Association Press, 2006,1–8


Nesme M, Kry P G, Jeřábková L, Faure F. Preserving topology and elasticity for embedded deformable models. ACM Transactions on Graphics, 2009, 28(3): 1–9 DOI:10.1145/1531326.1531358


Bickel B, Bächer M, Otaduy M A, Matusik W, Pfister H, Gross M. Capture and modeling of non-linear heterogeneous soft tissue. ACM Transactions on Graphics, 2009, 28(3): 1–9 DOI:10.1145/1531326.1531395


Faure F, Gilles B, Bousquet G, Pai D K. Sparse meshless models of complex deformable solids. In: ACM SIGGRAPH 2011 papers on-SIGGRAPH '11. Vancouver, British Columbia, Canada, New York, ACM Press, 2011, 30(4):76–79 DOI:10.1145/1964921.1964968


Chen D S, Levin D I W, Sueda S, Matusik W. Data-driven finite elements for geometry and material design. ACM Transactions on Graphics, 2015, 34(4): 1–10 DOI:10.1145/2766889


Chen D S, Levin D I W, Matusik W, Kaufman D M. Dynamics-aware numerical coarsening for fabrication design. ACM Transactions on Graphics, 2017, 36(4): 1–15 DOI:10.1145/3072959.3073669


Chen J, Bao H J, Wang T, Desbrun M, Huang J. Numerical coarsening using discontinuous shape functions. ACM Transactions on Graphics, 2018, 37(4): 1–12 DOI:10.1145/3197517.3201386


Chen J, Budninskiy M, Owhadi H, Bao H J, Huang J, Desbrun M. Material-adapted refinable basis functions for elasticity simulation. ACM Transactions on Graphics, 2019, 38(6): 161 DOI:10.1145/3355089.3356567


Johnson S G. The NLopt nonlinear-optimization package.


Stuart D A, Levine J A, Jones B, Bargteil A W. Automatic construction of coarse, high-quality tetrahedralizations that enclose and approximate surfaces for animation. In: Proceedings of Motion on Games-MIG '13. Dublin 2, Ireland, New York, ACM Press, 2019 DOI:10.1145/2522628.2522648


Barbič J, Sin F S, Schroeder D. VegaFEM Library,, 2012