CN101216950A - An elastic deformation simulation method based on linear differential operator - Google Patents

Abstract

The invention discloses an elastic deformation simulation method based on a linear differential operator. The invention comprises the following steps of using elastic energy based on the linear differential operator defined on the geometric attribute of a deformation gridding model to replace the traditional nonlinear energy based on continuum mechanics, carrying out modeling computation using Euler-Lagrange equation to compute the deformation process of the elastic object in the computer. Meanwhile, according to the calculation model, the invention optimizes the calculation process via further using the techniques of dimension reduction solution, airspace adaptive acceleration, etc, which allows the modeling process to reach the mutual or even real time efficiency and can simulate the elasticity of different materials and the deformation process of plastic objects. The invention reduces the calculation complexity, improves the calculation efficiency of the modeling solution, and keeps the living simulating effect via the model of transform solution based on the deformation energy form of the linear differential operator. The invention solves the rapid simulating problem in the deformation process of the elastic object in the real time or in a mutual virtual simulation system.

Description

A kind of elastic deformation simulation method based on the linear differential operator

Technical field

The present invention relates to the computer animation analogue simulation of elastic body object deformation, relate in particular to a kind of elastic deformation simulation method based on the linear differential operator.

Technical background

A lot of methods are often by simulating the result who obtains three-dimensional model deformation based on the elastic body of physics.Benefit is not only because can obtain the deformation result of the height sense of reality, and can obtain motion sequence to satisfy the needs of computer animation.The speed of simulation, stability, the error in the time of large deformation are wherein the most key several problems.Generally speaking, various algorithms are all in the balance of seeking between performance and the effect.And roughly can be divided into the method for linear elasticity model and the method for nonlinear elastic model based on the method for physical law.Choosing with the solution strategies of the equation of motion of elastic model is closely related, interacts.

In the physical simulation, the common method for expressing of deformable body is comprised the mass spring model, finite element model, no grid is expressed or the like.These methods all belong to Lagrangian type, that is to say that the sample record to object Ω (is called the material coordinate) in the object local coordinate system, and the motion of object has just been represented in the variation that the sampled point world coordinates takes place in simulation process.The mass spring model is the simplest a kind of.Sampled point characterizes the qualitative attribute of object, and resilient property is then expressed by the no quality spring that closes between the sampled point.Mostly very directly perceived simple based on the method for this thought of mass spring model, realize easily, and calculation cost is lower, but precision is relatively poor.

More accurate Finite Element Method based on continuum Model then is not simply object to be opened with particle is discrete, but uses tetrahedron, and elementary cells such as hexahedron are filled out the original.The attribute of object is in certain weighting (weighting function is commonly referred to type function, perhaps basis function B) for place, summit, unit thingness of unit internal representations.No grid expression has been removed the description of elementary cell topology in the finite element model, only pass through the discrete sampling point of object, and the type function on the sampled point is similar to the original.Though but the method precision that is based on this kind model is higher, calculate very complexity, extremely consuming time, can not be used for real-time or mutual analogue simulation and use.

Summary of the invention

The objective of the invention is at the deficiencies in the prior art, propose a kind of elastic deformation simulation method based on the linear differential operator.This method is used a kind of based on the elastic energy that is defined in the linear differential operator on the deformation grid model geometric attribute, substitute traditional non-linear energy, utilize Euler-Lagrangian equation of motion in computing machine, to carry out analog computation for the deformation process of elastomeric objects based on continuum mechanics.Simultaneously, at this kind computation model, further utilize technology such as dimensionality reduction is found the solution, space domain self-adapted acceleration that computation process is optimized, make simulation process can reach mutual even real-time efficient, and can simulate the elasticity of unlike material and the deformation process of plastic objects.

For realizing above-mentioned purpose, the technical solution used in the present invention comprises following four steps:

1) be defined in a class linear differential operator on the three-dimensional model geometric attribute, be used for describing three-dimensional model surface or inner distortion situation, different differentiating operators is corresponding to different simulate effects, and defines corresponding deformation energy at every kind of differentiating operator;

2) for above-mentioned differentiating operator, at first use the least square technology to solve the linear transformation of an optimum at each place, summit, use the technology of space domain self-adapted acceleration again, utilize the utmost point of bond quality weighting to decompose (polardecomposition) method, calculate the local rotation amount at each place, summit of three-dimensional model;

3) the local rotation amount that obtains deformation energy and the integrating step 2 that defines in the use step 1)), the process differentiate calculates stiffness matrix and the gradient matrix in the elastic analogy;

4) with step 1), 2), 3) in each component of asking for bring Euler-lagrange equation of motion into, use the dimensionality reduction solution technique, use implicit Euler method to carry out analog approach.

The described linear differential operational form that is defined on the three-dimensional model geometric attribute comprises Laplce (Laplacian) operator based on the summit, is used for simulating the motion deformation effect of entity; Based on the direction gradient operator on limit, be used for simulating shell deformation effect.

Described Laplace operator based on the summit, it is defined as follows:

$L_{i}x={\mathrm{\Σ}}_{j\∈N_i}{w}_{\mathrm{ij}}(x_j-x_i)$

L wherein _iFor being defined in i the differentiating operator on the summit, N _iThe all of its neighbor summit of expression node i.w _IjThe Laplace operator coefficient of the cotangent form on the expression undeformed model meshes, x _iBe the apex coordinate before i the summit undeformed, x is a unknown number, the vertex position after the expression deformation, and corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\left|\right|L_{i}x-R_i\left(x\right)d_i\left|\right|}^2$

Wherein V (x) is for being defined in the deformation energy on the whole model meshes, and K represents the topological connection relation of this grid, R _iExpression node i place is with respect to the part rotation of undeformed state, and scalar lambda is equivalent to Young modulus, describes the soft or hard degree of material, d _iThe Laplce coordinate of expression node i on the undeformed model meshes.

Described direction gradient operator based on the limit, it is defined as follows:

L _ijx＝w _ij(x _j-x _i)

L wherein _IjFor being defined in the differentiating operator on i summit and j the limit that the summit constituted, w _IjFor the limit (value of weight is relevant with the sampling density of model for i, weight j),

Corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\mathrm{\Σ}}_{(i,j)\∈K}{\left|\right|L_{\mathrm{ij}}x-R_i\left(x\right)d_{\mathrm{ij}}\left|\right|}^2$

D wherein _IjExpression limit (i, j) the Laplce's coordinate on the undeformed model meshes.

Described local rotation amount for each place, summit of three-dimensional model, at first use the least square technology to solve the linear transformation of an optimum at each place, summit, the utmost point decomposition method that re-uses the bond quality weighting is found the solution, and uses the technology of space domain self-adapted acceleration to be optimized to this process.

Described use least square technology solves the linear transformation of an optimum at each place, summit, utilize the linear transformation definition that problem is converted into the form of a least square energy of optimization and finds the solution and obtain analytical form as a result, thereby further use utmost point decomposition method to extract rotational invariants.

The technology of the space domain self-adapted acceleration of described use is adaptively chosen representational node, calculates its local rotation, then it is diffused on other summits, and algorithm carries out breadth First traversal in the mode of greediness to node from the minority seed points.In ergodic process, automatically identify representation node.

Described breadth First traversal, in the ergodic process, need to judge whether a rotation matrix can represent the part rotation of other nodes, in the transformation matrix process of comparison node, adopt " rotation is sensitive " criterion, ignore flexible, cut mistake, only rotating part has relatively improved efficient and accuracy relatively greatly.

Described application dimensionality reduction solution technique is simplified computation model, and has transformed the form of Euler-lagrange formula, and the variable dimension that needs are found the solution reduces greatly, re-uses implicit Euler method and carries out analog approach.

Described dimensionality reduction solution technique, too high at elastic body simulation degree of freedom, cause problem scale excessive and find the solution the not high problem of stability, utilization dimensionality reduction strategy projects to the higher-dimension variable that needs to find the solution in the deformation simulation in the linear space of low-dimensional.Comprising two kinds of methods: based on the method for the distortion sample of input, for given sample, use the mass-PCA choice of technology go out most important base vector as dimensionality reduction after the orthogonal basis in deformation space; Based on method of mode analysis, for the occasion that does not have the deformation sample,, the deformation space is decomposed by finding the solution generalized eigenvalue problem, choose little eigenwert characteristic of correspondence vector and be the base in deformation space.

The present invention compares with background technology, and the beneficial effect that has is:

The proposition of this method novelty based on the deformation energy form of linear differential operator, by the model that conversion is found the solution, reduced the complexity of calculating, improved the counting yield of analog approach greatly, kept simulate effect very true to nature simultaneously.This method solved in real time or in the mutual level dummy emulation system for the quick problem of modelling of the deformation process of elastomeric objects, can be used to numerous areas such as analog simulation, education demonstration, Entertainment.

Description of drawings

Fig. 1 is a process flow diagram of the present invention.

Fig. 2 is the surface mesh simulate effect figure of different differentiating operators.

Fig. 3 is the process flow diagram of the greedy method of diffusion of the breadth First that adopted in the space domain self-adapted speed technology.

Fig. 4 is the distribution plan of self-adaptation node in the space domain self-adapted speed technology.

Fig. 5 is based on the effect synoptic diagram of the analysis dimension reduction method of sample in the dimensionality reduction solution technique.

Fig. 6 is based on the effect synoptic diagram of modal analysis method in the dimensionality reduction solution technique.

Embodiment

The invention will be further described below in conjunction with drawings and Examples.

As shown in Figure 1, a kind of elastic deformation simulation method based on the linear differential operator that the present invention proposes comprises that the calculating of the definition of differentiating operator, extraction, stiffness matrix and gradient matrix that local rotation amount is located on the summit and the equation of motion four steps such as find the solution.

Now specifically introduce four steps of this method:

1) definition of differentiating operator

The definition of the present invention on the three-dimensional model geometric attribute one class linear differential operator, be used for describing three-dimensional model surface or inner distortion situation, different differentiating operators is corresponding to different simulate effects, and defined corresponding deformation energy at every kind of differentiating operator.Comprising: Laplce (Laplacian) operator based on the summit is used for simulating the motion deformation effect of entity; Based on the direction gradient operator on limit, be used for simulating shell deformation effect.Fig. 2 has provided the surface mesh simulate effect of different operators.Fig. 2 a is a master mould; Fig. 2 b is based on the deformation effect of the Laplacian operator on summit; Fig. 2 c. is based on the deformation effect of direction gradient operator.

Based on the Laplacian operator on summit, it is defined as follows:

$L_{i}x={\mathrm{\Σ}}_{j\∈N_i}{w}_{\mathrm{ij}}(x_j-x_i)$

L wherein _iFor being defined in i the differentiating operator on the summit, N _iThe all of its neighbor summit of expression node i.w _IjThe Laplacian operator coefficient of the cotangent form on the expression undeformed model meshes.x _iIt is the apex coordinate before i the summit undeformed.X is a unknown number, the vertex position after the expression deformation.

Corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\left|\right|L_{i}x-R_i\left(x\right)d_i\left|\right|}^2$

Wherein V (x) is for being defined in the deformation energy on the whole model meshes, and K represents the topological connection relation of this grid.R _iExpression node i place is with respect to the part rotation of undeformed state.Scalar lambda is equivalent to Young modulus, describes the soft or hard degree of material.d _iThe Laplacian coordinate of expression node i on the undeformed model meshes.

Based on the direction gradient operator on limit, it is defined as follows:

L _ijx＝w _ij(x _j-x _i)

L wherein _IjFor being defined in the differentiating operator on i summit and j the limit that the summit constituted.w _IjBe limit (i, weight j).The value of weight is relevant with the sampling density of model.

Corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\mathrm{\Σ}}_{(i,j)\∈K}{\left|\right|L_{\mathrm{ij}}x-R_i\left(x\right)d_{\mathrm{ij}}\left|\right|}^2$

D wherein _IjExpression limit (i, j) the Laplacian coordinate on the undeformed model meshes.

2) extraction of local rotation amount

In order to guarantee the rotational invariance of deformation process, must calculate each place, summit distortion local rotation amount afterwards, thereby could correctly calculate deformation energy.For above-mentioned differentiating operator, the present invention at first uses the least square technology to solve the linear transformation of an optimum at each place, summit, use the technology of space domain self-adapted acceleration again, utilize the utmost point of bond quality weighting to decompose (polar decomposition) method, calculate the local rotation amount at each place, summit of three-dimensional model.

Utilize the linear transformation definition that the form that problem is converted into a least square energy of optimization is as follows:

Make r _iOriginal position during for the node i undeformed, x _iBe deformation posterior nodal point location definition:

r _ij＝m _j(r _j-r _i)

x _ij＝m _j(x _j-x _i)

Then at node i place minimization ∑ _{J ∈ N (i)}‖ x _Ij-A _ir _Ij‖ ²Obtain the optimum linear conversion $A_i=A_i^{\mathrm{xr}}{A}_i^{\mathrm{rr}}$ , wherein:

$A_i^{\mathrm{xr}}=\left(\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{x}_{\mathrm{ij}}{r}_{\mathrm{ij}}^T\right)$ $A_i^{\mathrm{rr}}=\left(\underset{j\∈N\left(i\right)}{\mathrm{\Σ}}{r}_{\mathrm{ij}}{r}_{\mathrm{ij}}^T\right)$

We just can use the utmost point decomposition method of bond quality weighting to this linear transformation A then _iThereby decompose the local rotation amount that obtains each place, summit.

But it is very big all to carry out such operational computations expense for each summit.The present invention utilizes the technology of space domain self-adapted acceleration that computation process is carried out following improvement: adaptively choose representational node, calculate its local rotation, then it is diffused on other summits.Algorithm carries out breadth First traversal in the mode of greediness to node from the minority seed points.In ergodic process, automatically identify representation node.Visible Fig. 3 of idiographic flow of the strategy of the greedy diffusion of this breadth First.In above-mentioned diffusion process, of paramount importance step is to judge that can a matrix R accurately substitute the part rotation of other nodes.The present invention adopts " rotation is sensitive " criterion, ignore flexible, cut mistake, only rotating part has relatively improved efficient and accuracy relatively greatly.This judgment criterion is as follows:

Calculate ‖ R _i ^TA _j ^Xr-(R _i ^TA _j ^Xr) ^T‖ ²If the value of this formula is less than certain threshold value, then rotation R _iBe diffused into node j from node i.Shown among Fig. 4 that adaptive representation node distributes and deformation effect figure.

3) calculating of stiffness matrix and gradient matrix

At first the differentiating operator matrix is reset, it is as follows to be organized into more compact form:

$L=\left(\begin{array}{c}{L}_1\\ L_2\\ \·\\ \·\\ \·\end{array}\right)$ $R\left(x\right)=\left(\begin{array}{ccccc}{R}_1\left(x\right)& & & & \\ & R_2\left(x\right)& & & \\ & & .& & \\ & & & .& \\ & & & & .\end{array}\right)$ $d=\left(\begin{array}{c}{d}_1\\ d_2\\ \·\\ \·\\ \·\end{array}\right)$

Wherein L is the matrix of expression linear differential operator, and the unique piecemeal diagonal matrix of R (x), its each piecemeal are represented the local rotational component on the corresponding vertex, and d represents the vector that all micro components are lined up.

Then the deformation energy The Representation Equation is as follows:

$V\left(x\right)=\frac{1}{2}{\left|\right|\mathrm{Lx}-R\left(x\right)d\left|\right|}^2$

This deformation energy is calculated its gradient and black plug (Hesse) matrix obtains:

$\frac{\∂}{\∂x}V\left(x\right)=\mathrm{\λ}{(L-\frac{\∂R\left(x\right)}{\∂x}d)}^T(\mathrm{Lx}-R\left(x\right)d)\≈{\mathrm{\λL}}^T(\mathrm{Lx}-R\left(x\right)d)$

$\frac{{\∂}^2}{\∂x^2}V\left(x\right)=\mathrm{\λ}{L}^T(L-\frac{\∂R\left(x\right)}{\∂x}d)\≈\mathrm{\λ}{L}^{T}L$

Thus, obtain the approximate quantity λ L of stiffness matrix ^TL.

4) equation of motion is found the solution

Use the dimensionality reduction solution technique, transform the form of Euler-lagrange formula.Bring into finding the solution each component that obtains in the above process again, use implicit Euler method to carry out analog approach.Equation form is as follows:

${\mathrm{\Φ}}^T(M+\mathrm{hD}+h^{2}K)\mathrm{\Φ\Δ}\stackrel{\·}{z}=h{\mathrm{\Φ}}^T(f_{\mathrm{ext}}-D\stackrel{\·}{x}-\mathrm{hK}\stackrel{\·}{x}-\frac{\∂}{\∂x}V\left(x\right))$

Wherein Ф represents the orthonormal basis matrix in deformation space, and z is the parameter of deformation in the dimensionality reduction space.The present invention adopts two kinds of methods to obtain the base in deformation space.A kind of method that is based on the distortion sample of input: for given m sample { s _j} _J=1 ^m, at first to deformation vector { s _j-r} _J=1 ^mUse the mass-PCA technology, select the orthogonal basis of most important base vector then as dimensionality reduction deformation space.Provided two groups of users' given sample among Fig. 5, and the deformation result who obtains by these samples.

The another kind of method that obtains dimensionality reduction orthogonal space base is based on method of mode analysis, is applicable to the occasion that does not have the deformation sample.At first find the solution generalized eigenvalue problem M Ψ Λ=K Ψ, each of matrix Ψ is classified the generalized eigenvalue vector as, and the diagonal entry of diagonal matrix Λ is an eigenwert.Choose little eigenwert characteristic of correspondence vector and be the base in deformation space, they have represented the part that accounts for main energy of elastic movement medium and low frequency.Fig. 6 represents the synoptic diagram based on modal analysis method, and the left side is 8 non-trivial bases in the deformation space that obtains of model analysis, and the right is the deformation effect.

Claims (10)

1. elastic deformation simulation method based on the linear differential operator is characterized in that comprising following four steps:

1) be defined in a class linear differential operator on the three-dimensional model geometric attribute, be used for describing three-dimensional model surface or inner distortion situation, different differentiating operators is corresponding to different simulate effects, and defines corresponding deformation energy at every kind of differentiating operator;

2) for above-mentioned differentiating operator, at first use the least square technology to solve the linear transformation of an optimum at each place, summit, use the technology of space domain self-adapted acceleration again, utilize the utmost point decomposition method of bond quality weighting, calculate the local rotation amount at each place, summit of three-dimensional model;

3) the local rotation amount that obtains deformation energy and the integrating step 2 that defines in the use step 1)), the process differentiate calculates stiffness matrix and the gradient matrix in the elastic analogy;

4) with step 1), 2), 3) in each component of asking for bring Euler-lagrange equation of motion into, use the dimensionality reduction solution technique, use implicit Euler method to carry out analog approach.

2. a kind of elastic deformation simulation method according to claim 1 based on the linear differential operator, it is characterized in that: the described linear differential operational form that is defined on the three-dimensional model geometric attribute, comprise Laplace operator, be used for simulating the motion deformation effect of entity based on the summit; Or, be used for simulating shell deformation effect based on the direction gradient operator on limit.

3. a kind of elastic deformation simulation method based on the linear differential operator according to claim 2 is characterized in that: described Laplace operator based on the summit, and it is defined as follows:

$L_{i}x={\mathrm{\Σ}}_{j\∈N_i}{w}_{\mathrm{ij}}(x_j-x_i)$

L wherein _iFor being defined in i the differentiating operator on the summit, N _iThe all of its neighbor summit of expression node i.w _IjThe Laplace operator coefficient of the cotangent form on the expression undeformed model meshes, x _iBe the apex coordinate before i the summit undeformed, x is a unknown number, the vertex position after the expression deformation, and corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\left|\right|L_{i}x-R_i\left(x\right)d_i\left|\right|}^2$

Wherein V (x) is for being defined in the deformation energy on the whole model meshes, and K represents the topological connection relation of this grid, R _iExpression node i place is with respect to the part rotation of undeformed state, and scalar lambda is equivalent to Young modulus, describes the soft or hard degree of material, the Laplce coordinate of expression node i on the undeformed model meshes.

4. a kind of elastic deformation simulation method based on the linear differential operator according to claim 2 is characterized in that: described direction gradient operator based on the limit, and it is defined as follows:

L _ijx＝w _ij(x _j-x _i)

L wherein _IjFor being defined in the differentiating operator on i summit and j the limit that the summit constituted, w _Ij(value of weight is relevant with the sampling density of model for i, weight j) for the limit;

Corresponding deformation energy definition form is as follows:

$V\left(x\right)=\frac{\mathrm{\λ}}{2}{\mathrm{\Σ}}_{i\∈K}{\mathrm{\Σ}}_{(i,j)\∈K}{\left|\right|L_{\mathrm{ij}}x-R_i\left(x\right)d_{\mathrm{ij}}\left|\right|}^2$

D wherein _IjExpression limit (i, j) the Laplce's coordinate on the undeformed model meshes.

5. a kind of elastic deformation simulation method according to claim 1 based on the linear differential operator, it is characterized in that: described local rotation amount for each place, summit of three-dimensional model, at first use the least square technology to solve the linear transformation of an optimum at each place, summit, the utmost point decomposition method that re-uses the bond quality weighting is found the solution, and uses the technology of space domain self-adapted acceleration to be optimized to this process.

6. a kind of elastic deformation simulation method according to claim 5 based on the linear differential operator, it is characterized in that: described use least square technology solves the linear transformation of an optimum at each place, summit, utilize the linear transformation definition that problem is converted into the form of a least square energy of optimization and finds the solution the analytical form that obtains the result, thereby further use utmost point decomposition method to extract rotational invariants.

7. a kind of elastic deformation simulation method according to claim 5 based on the linear differential operator, it is characterized in that: the technology of the space domain self-adapted acceleration of described use, adaptively choose representational node, calculate its local rotation, then it is diffused on other summits, algorithm carries out breadth First traversal in the mode of greediness to node from the minority seed points, in ergodic process, automatically identify representation node.

8. a kind of elastic deformation simulation method according to claim 7 based on the linear differential operator, it is characterized in that: described breadth First traversal, in the ergodic process, need to judge whether a rotation matrix can represent the part rotation of other nodes, in the transformation matrix process of comparison node, adopt " rotation is sensitive " criterion, ignore flexible, cut mistake, only compare rotating part, improved efficient and accuracy relatively greatly.

9. a kind of elastic deformation simulation method according to claim 1 based on the linear differential operator, it is characterized in that: described analog approach, use the dimensionality reduction solution technique, simplify computation model, and transformed the form of Euler-lagrange formula, the variable dimension that needs are found the solution reduces greatly, re-uses implicit Euler method and carries out analog approach.

10. a kind of elastic deformation simulation method according to claim 9 based on the linear differential operator, it is characterized in that: described dimensionality reduction solution technique, too high at elastic body simulation degree of freedom, cause problem scale excessive and find the solution the not high problem of stability, utilization dimensionality reduction strategy projects to the higher-dimension variable that needs to find the solution in the deformation simulation in the linear space of low-dimensional.Comprising two kinds of methods: based on the method for the distortion sample of input, for given sample, use the mass-PCA choice of technology go out most important base vector as dimensionality reduction after the orthogonal basis in deformation space; Based on method of mode analysis, for the occasion that does not have the deformation sample,, the deformation space is decomposed by finding the solution generalized eigenvalue problem, choose little eigenwert characteristic of correspondence vector and be the base in deformation space.

Images