1. Introduction
Symmetry is a ubiquitous phenomenon in nature and a very essential rule for artistic design. Mathematically, it encodes some invariance of given structures under certain transformations. In practice, symmetry reflects both aesthetics considerations and functionality purposes. Because of its importance, symmetry has been widely used in many geometry processing procedures, such as remeshing, editing, repairing, etc, as well as the virtual reproduction of animations with symmetry for visual realism and artistic harmony.
However, at the stage of simulation, the symmetry may be lost. Some models used in computer animation come from 3D scanning or other geometry modeling software. On the one hand, limited by data acquiring techniques and real-world noises, the symmetry of a shape may be broken when transformed into digital representations. More importantly, even if the shape is in perfect symmetry, e.g., modeled by software, the mesh tessellated from the shape could be asymmetric. The elastic simulation based on the asymmetric mesh may result in asymmetric behavior. In other words, the displacement or acceleration may be not symmetric even if the force field is in symmetry (see
Figure 1). To reduce the asymmetric artifacts, increasing the mesh density may help, but the additional computational cost is not affordable in some real-time applications (e.g., game, virtual reality). If the mesh has some asymmetrically distributed badly shaped elements, which brings strong shear locking, we also notice that refining the mesh by subdivision cannot easily reduce the asymmetrical behavior (see
Figure 1).
To get rid of such undesired asymmetry artifacts, a natural idea is to generate a symmetric mesh, and a few works (e.g., [
1]) do make such attempts. However, they are not generally applied and cannot be used for generating the desired volumetric meshes. It is worth noticing that it is even impossible to generate meshes with the desired symmetry sometimes. For instance, no structural hexahedral mesh can have the correct symmetry at the rotational axis of a fan with five symmetric blades.
There are other possible remedies. One can impose constraints or use post-processing to modulate the elastic simulation so that the displacement or even acceleration is symmetric. However, such a way introduces artificial ingredients into the simulation, and the constraints or the post-processing are not easy to design for general situations. Roughly speaking, such a way is very restricted and not applicable in general elastic simulation involving asymmetric boundary conditions. A more flexible method is to refine the mesh or use a high order basis function to reduce the impact of asymmetry in the mesh. This strategy results in a larger or denser system matrix, and the additional cost may not be affordable for some applications that require intensive user interactions and rapid feedbacks.
Considering the above-mentioned issues, we instead propose to optimize the material to compensate the error coming from asymmetric meshes. Unlike refinement approaches, our method keeps the problem size and matrix sparsity. Consequently, numerical performance will not be notably affected while desired symmetry behaviors can be nicely preserved with the optimized materials. Moreover, once the material gets adjusted, it can be generalized to solve a series of problems defined over this object without any necessities for specifically designed constraints or post-processing, thus saving enormous cost on computational or operational efforts and allowing for better interactivity for editing or design tasks.
Technically, there are two challenges to be addressed. Firstly, a clear mathematical definition of symmetric behavior is required. Secondly, a reasonable symmetry-deviation measurement is necessary—because even if the material gets optimized, it is generally not possible to have symmetric behavior in any case. To address these two challenges, we make the following contributions:
Both static and kinematic symmetry behaviors are considered. Under a symmetric force field, we ask the displacement of the static equilibrium status and acceleration in kinematic simulation to be symmetry.
We find that all the above symmetric vector fields fall in the kernel of a linear operator associated with the symmetry. Therefore, the asymmetric error can be efficiently measured via a simple norm.
Taking the common assumption used in modal analysis that low frequency modes dominant the behavior, we formulate the physical equation in the above kernel space and take its leading eigenvectors as the training force fields so that the resulting material can be generally applied.
To justify the proposed method, we show results on various meshes with rotational and mirror shape symmetries, as well as an illustration on symmetric behavior in dynamic elastic simulation. The asymmetric error generally deceases greatly after using an optimized material, which clearly shows the benefits of our method.
4. Discretization and Numerical Method
In this section, we show our implementation for the above formulations. We show how to apply the aforementioned method to an model discretized by a tetrahedral mesh composed of a node set and an element set with piece-wise linear shape functions. Other discretization strategies, such as higher order FEM or even meshless methods, are also feasible.
To discretize the key equation Equation (
7), some points
in
are sampled, and the integration is approximated by the weighted sum of
. To be specific, the sampling points are the boundary nodes of
in
, i.e.,
. The weight for each sampling point is the area of the Voronoi cell in
. It should be noticed that we indeed just evaluate the integration on the boundary of the mesh, but it works well in our experiments. The set
contains all the points
symmetric to
. If
happens to fall outside of
, we use its nearest point on
to replace it.
Denoting the flattened vectors of all nodal coordinates in
, and
as
,
, and
, respectively, we can have two sparse matrices
to interpolate the positions of the points, i.e.,
and
. The symmetry of these two point sets can be verified by:
where
and
. The linear operator
in Equation (
7) is then discretized as a
matrix:
Following the common linear FEM strategy, a 3D vector field
z, which can be a displacement field
u, acceleration field
a, or a force field
f, is encoded on each node and linearly interpolated in each element, and thus can be represented as
. The vectors on the above points can be interpolated by
. As a consequence, the integration in Equation (
7) can be approximated as:
where
is a
diagonal matrix whose entries are the aforementioned weight for each sampling point.
The discretization of the elastic model and other parts is in the common way. The scale of Young’s modulus
s is piece-wise constant in each element, and thus can be represented as a vector
. The global stiffness matrix
is assembled from local stiffness matrices computed on each tetrahedral element. Because scaling Young’s modulus of the original material is equivalent to scaling the local stiffness matrix, we have:
where
is the scale in the
e-th element. The stiffness matrix
for element
e is pre-computed with the input Young’s modulus. The scale of mass is encoded on each node, thus represented as a vector
as well. We use the lumped mass matrix, and the entries related to the
n-th node in this diagonal matrix
have the value
where
stands for the input mass of this node. The energy to fix symmetric points is simply reduced to:
where
is a
sparse matrix to interpolate the displacement of
and
from
. The stiffness matrix for these constraints is
.
Putting all of them together, Equation (
11) becomes a linear system in
dimensions:
The training force field
is in the kernel of
. Moreover, those forces of the nodes that are not involved by
are all
, which intuitively means that we do not apply any force on such nodes not in
and
or algebraically meaning that if the column of
is all zeros, then the corresponding component of
is zero. Actually, we are more interested in some symmetric forces, which implies that this force is in the kernel space of
S. A set of bases for
can be extracted by singular value decomposition of
. Given a matrix
with columns composed of such bases, we solve the following general eigenvalue problem for
associated with non-trivial eigenvalues:
is similar to
, but composed of the Voronoi volume for every node in
. The training force fields simply become
, where
are the columns of
. For better understanding, We visualize the shape of the training forces of different models in
Section 5.3. The weight of asymmetric measurement for
is
, and the number of forces
N is set to 40 experimentally in all of our experiments. One can set
as
or
as
, but this will significantly increase the computational cost of the singular value decomposition here.
The asymmetric energy term
F in Equation (
17) is thus discretized as:
and the regularization term is:
where
is an
diagonal matrix with entries
being the volume of the
i-th element in mesh
, and
here is a constant
-D vector whose entries are all ones.
Finally, the optimization problem becomes:
Note that when we are further attempting to symmetrize kinematic behavior, the objective function can be thus replaced by:
in conjunction with the additional regularization term and box constraints on variable
. There is no need to alternatively optimize the stiffness scaling
and mass scaling
, because they are indeed decoupled and used to symmetrize the static equilibrium and kinematic behavior, respectively. The stiffness scaling can be simply fixed after getting optimized in the equilibrium step for efficiency, which is also helpful to keep an invariant static equilibrium state as the object finally stops moving after a long time.
In addition, some user-controlled parameters are summarized here: , and h are set to , , 1, and , respectively.
Numerical Method
Because the material to compensate the asymmetric mesh is not smoothly changed in the general situation, existing acceleration techniques like dimension reduction [
14] are not applicable. Instead, we solve the material optimization problem in full space.
The Hessian of the objective function in Equation (
28) is a dense matrix, which leads to very high computational cost by the common Newton method. Therefore, we take the L-BFGS method provided by the ALGLIB library [
27] with box constraints because it just needs the gradient of the objective function derived in
Appendix A. During the optimization, we set
as the initial guess and take a simple diagonal preconditioner provided in
Appendix B to accelerate the convergence. We use multiple conditions for termination: if one of the following conditions satisfies, i.e., the change of the objective value,
or the change of the step value,
, or the norm of the gradient,
is less than the corresponding threshold, or the number of iterations reaches the maximum, the iteration stops. In practice, The first three thresholds
,
, and
are set to
, and the maximal iteration number is set to 500 in all of our experiments.
5. Results
Unless otherwise specified, we use default parameters mentioned in the last section in the optimization process for all of the following experiments.
5.1. Optimized Symmetric Deformation
In this section, we present side-by-side comparisons of the deformations before and after material optimization on models with reflectional or rotational symmetry.
In
Figure 3, for a cuboid model, we take points
on its reflection plane as fixed and apply the same forces on its two ends (as illustrated in the first row). In the first column, the deformations before and after optimization are presented when the same downward forces are applied on both ends. It is clear to see that the deformation result after optimization becomes more symmetric. The effectiveness of the optimization is also shown in the second row when the outward forces are imposed. Furthermore, the histogram in
Figure 3 for the distribution of the symmetry error on all the nodes clearly demonstrates the effectiveness.
Similar results can be obtained for models with rotational symmetry. The cross model in
Figure 4 has four symmetric parts. Just like the fan model in
Figure 2, it has three multiple
rotational symmetries. We fix points on the innermost layer (as shown by the red points in the first row) and illustrate its deformation results under two symmetric forces before/after optimization. As is shown in the figure, the deformations under the downward or the spin force after optimization are much more symmetric than those before optimization as well.
It is worth noticing that although being simple in geometry, these two examples clearly show the benefits of our material adaption regarding the deformation symmetry. For ever more complex models, things could get even worse, where neither discrete shapes nor the tessellations can easily stay perfectly symmetric in the discretization. A rough, irregularly sampled, or even non-manifold input surface mesh could seriously destroy the regularity of the resulting volumetric meshes, further making the solutions on the generated meshes look hardly close to symmetric, while our method notably overcomes existing defects in meshing.
Furthermore, to conduct quantitative assessment on the results, we also show the symmetric displacement error on each sampling point in one of the symmetric parts (e.g.,
in Equation (
7)). The symmetric displacement error
on a sampling point
is defined as:
In
Figure 5, we compare both the deformation results and the distribution of
on a few models before/after optimization. In this figure, all symmetric forces are applied on the ends of models with the center part of the models being fixed. In such a setting, sampling points near the end of the model generally have larger symmetric displacement error. As for the symmetric displacement error
, it tends to be smaller after optimization, which leads to a better deformation result (more symmetric).
5.2. Optimized Material Distribution
In this section, we look into the optimized material. The optimized material can be represented by the stiffness scale
on each element. In
Figure 6, we visualize the scale field
on the cuboid model (
Figure 3) and the cross model (
Figure 4). Note that these two models have very different tessellations in the symmetric parts. It can be observed that, in the symmetric parts, big elements tend to have a small stiffness scale, which could compensate the numerical stiffness. On the contrary, small elements tend to have a large stiffness scale, making them more stiff. Overall, the optimized material makes the deformation result symmetric under symmetric forces, as was illustrated in the previous section.
Most of our previous experiments fix points near the reflection plane or rotation axis of a model, which is common in practice. However, our method does apply to other scenarios as well. For example, for the cuboid model, we can also fix points on its ends and apply symmetric forces to the other part of the model. As shown in
Figure 7, the sampling points on the end of the cuboid (red sphere in (a)) are fixed, while downward force is applied to the whole model. In this case, the deformation result after material optimization also becomes more symmetric. From the comparison results of
Figure 6 and
Figure 7, we can see that different fixed points lead to a different optimized material. This implies that the material optimized for one boundary condition may not be suitable for others. As a limitation, this will be discussed at the end of this paper.
5.3. Kernel Forces Visualization
To make it easy and intuitive to understand the force at different modes, we visualize some of these modes in
Figure 8.
5.4. Performance
The performance of our algorithm is investigated in this section. In all of our experiments, the optimization energy in Equation (
28) drops quickly in the first few iterations and mostly converges in less than 350 iterations, as is shown in
Figure 9.
We also count and analyze the time consumption of our algorithm. In our implementation, there are two time-consuming stages: one is the generation of the training force
, and the other is the numerical optimization (Equation (
28)). Suppose that one is able to know what kinds of symmetric forces are mostly used in practice, then the optimization can only focus on the forces of interests, which saves the time-consuming effect of computing many kernel vectors in the symmetry subspace. We show the statistic results in
Table 1. These results were evaluated on a PC with an Intel®Core™i5-6300HQ CPU @ 2.30 GHz × 4, and 8 G memory.
Although such a precomputation step is inefficient, extra cost can be saved in the following stages of numerical optimization. This reveals the fact that our method for symmetrizing solutions does not need to have sharply increasing cost due to mesh or basis refinement. Thus, it is particularly useful to process batches of tasks, e.g., generating a many-frame animation of nonlinear dynamics, as we will show in the next subsection.
5.5. Symmetrizing Dynamics
As mentioned, our idea to adjust the material is not only applicable to elastic stiffness for pursuing symmetric equilibrium statics, but is also able to symmetrize the dynamics by incorporating mass compensation. Once the mass and stiffness have been optimized, the physical model can be applied to simulate a sequence of symmetry-preserving animations without the need for frame-by-frame manual editing. In this subsection, we show that our method provides a preliminary method for the dynamics simulation. Indeed, our method just “symmetrizes” the dynamic behavior around the rest state. Integrating the co-rotational method, we can achieve some natural results under large deformation. However, general non-linear elastic models will bring big challenge to our method, which is reserved for future exploration.
Firstly, we will show how such an additional mass rescaling helps to symmetrize the dynamic simulations in terms of shapes. Although the material optimization in our method is conducted based on the rest-state stiffness, we found that the material scalings are also applicable to co-rotational elasticity, which has the same underlying linear strain-stress relationship, but properly gets rid of rotations when computing the strain from displacements. In
Figure 10, it is clear that directly simulating a bird flying with the original mesh causes obvious asymmetry along the dynamic sequence because of the asymmetry of the input mesh. Unlike static deformation, such discretization induced artifacts may become severe since they could be accumulated in the dynamic evolution. In contrast, our method managed to produce a far more symmetric co-rotational simulation in a completely automatic way.
Figure 11 shows that our method can preserve the symmetry in the dynamic trajectory. We apply the same load on each wing of a fan model for an evenly spaced time duration as the initial condition. As illustrated in
Figure 11, the dynamic trajectories show notable phase error without material optimization. The results are refined when the stiffness is optimized, and it keeps a significantly better symmetry in “phase” while both the mass and stiffness are optimized. Note that it is difficult or tedious for traditional methods to design special constraints or post-processing for this simulation.