CN104318599A - High-precision fluid animation modeling method based on geometrical features - Google Patents

Abstract

The invention discloses a high-precision fluid animation modeling method based on geometrical features. The high-precision fluid animation modeling method includes the following steps: initializing initial conditions and boundary conditions in a fluid physical model; obtaining the parameters which are required by calculation of the geometrical features and include an intermediate speed field, a fluid speed field in a current time step and a symbol distance field function; building a geometrical model of a fluid animation, and carrying out automation location, quantization and tracking on detail sensitive areas in the fluid animation; according to a feature function of the geometrical features, controlling particle sets of the physical model, and coupling the geometrical model with the physical model; extracting triangular grid surfaces representing the fluid surface and fluid sprays, and carrying out reality sense rendering on the triangular grid surfaces to generate the frame of fluid animation. By means of the high-precision fluid animation modeling method, the detailed effects of the high-precision fluid water surface, the high-precision fluid sprays and the like can be effectively achieved, automation generation of the high-precision fluid animation is achieved, and the high-precision fluid animation modeling method has the advantages of being high in precision and high in efficiency.

Description

A kind of high precision fluid animation modeling method based on geometric properties

Technical field

The present invention relates to a kind of high precision fluid animation modeling method based on geometric properties.

Background technology

Fluid animation has a wide range of applications in virtual reality, cartoon special effect, interactive game, architectural design, product introduction, medical science, military simulation emulation etc.In these fluid animations, visual vivid is the most basic pursuit.But a lot of fluid animation is all that animation teacher makes by hand in tool software at present.Because fluid motion is very complicated, even have the animation teacher of rich experiences, be also difficult to the animation that arbitrary design is wanted.In addition fluid animation manufacturing speed is very slow, cost is high, is typical " labor-intensive " industry.Therefore, how automatically generation has the problem that the fluid animation enriching details is a very challenging property.Traditional fluid animation method exists that calculation cost is high, efficiency is low and the problem that fluid details is easily lost, and is difficult to the high precision fluid animation generating visual vivid.

Summary of the invention

The present invention is in order to solve the problem, propose a kind of high precision fluid animation modeling method based on geometric properties, this method the method is by generating the high precision fluid animation robotization of visual vivid, its treatment effeciency in the applications such as computer graphics, computer animation, special effect making, Entertainment of effective raising, making the transition for promoting animation industry from " labor-intensive " to " technology-intensive type " further, expanding its application and scope based theoretical.

To achieve these goals, the present invention adopts following technical scheme:

Based on a high precision fluid animation modeling method for geometric properties, comprise the following steps:

(1) starting condition in convection cell physical model and boundary condition carry out initialization;

(2) from the physical model of fluid, obtain the fluid velocity field in the midrange speed field required for geometric properties calculating, current time step and symbolic distance field function parameter;

(3) build the geometric model of fluid animation, robotization location, quantification and tracking are carried out in the details sensitizing range in convection cell animation;

(4) according to the fundamental function of geometric properties, the particle collection of physical model is controlled, coupling geometric model and physical model;

(5) extract the triangular mesh surface representing flow surface and fluid spray, and Realistic Rendering is carried out to it, generate a frame fluid animation;

(6) check the making whether completing whole fluid animation, if do not completed, then return step (2), proceed to make.

In described step (1), the initialized step of physical model comprises:

(1-1) the fluid velocity field u in initial background grid;

(1-2) initialization represents the symbolic distance field function φ of flow surface;

(1-3) initialization convection cell surface carries out detecting the particle collection P corrected _cwith the particle collection P representing fluid spray _s;

(1-4) initialization border, namely uses marker bit to represent the type of each grid cell in fluid calculation territory; Need, by solid obstacle by computational fields grid cell size discrete, to be set to 1 by the grid cell marker bit that solid occupies after discretize, otherwise to be set to 0 for this reason.

In described step (2), the concrete grammar obtaining geometric properties from the physical model of fluid is as follows:

(2-1) according to the fluid velocity field u in initial background grid, particle collection P _sin speed calculate midrange speed field in background grid;

(2-2) according to midrange speed field, by solving the pressure Poisson equation of coupling Euler and Lagrange two kinds of models, the fluid velocity field u in current time step is obtained ⁿ⁺¹;

(2-3) according to obtained current time step fluid velocity field, the symbolic distance field function φ after upgrading is calculated ⁿ⁺¹.

In described step (3), set up the geometric model of fluid animation, the details sensitizing range in convection cell animation carry out robotization location, quantification and tracking basic step as follows:

(3-1) build envelope closed manifold, is formed by sealing fluid interface or nonocclusive fluid interface and computational fields border, seal closed manifold border and surrounding 3 layers of grid are called narrowband region;

(3-2) definition local axis (Local Medial Axis, LMA) is for being positioned at arrowband Ω in envelope closed manifold _bon the track in the inscribed circle center of circle;

(3-3) according to the symbolic distance field function φ obtained in Step2 ⁿ⁺¹, use divergence value determining method to calculate LMA;

(3-4) according to local axis, definition location quantization function (Location and Quantization Function, LF), for stream shape borderline net point, its functional value is given threshold value λ and the smaller value of this point to locally half-breadth value; In narrowband region, the location quantization function value of other net points is that this point adds the location quantization function value of shape border, upper reaches corresponding point to the bee-line flowing shape border; Location quantization function value for other net points on computational fields is set to a large constant, and location quantization function can realize smooth location to deep camber, elongated and edge details sensitizing range and quantification;

(3-5) according to the definition of location quantization function, by using the method for KD retrieval, the location quantization function value of Fluid Computation net point on the surface;

(3-6) according to the definition of location quantization function, by using the location quantization function value extrapolation of fast marching algorithms convection current shape frontier point, the location quantization function value of other net points in narrow-band-domain is calculated;

(3-7) according to the definition of location quantization function, the location quantization function value in setting background grid beyond narrow-band-domain is large constant, is appointed as N _max× Δ h, wherein N _maxfor computational fields is in the maximal value of x, y, z three axially meshes number, Δ h is the size of grid cell.

(3-8) deformation extent on accumulation Fluid Computation surface in each time step, if deformation extent is less than setting threshold value, then needs to carry out convection current renewal to location quantization function; Otherwise need to reinitialize location quantization function.

In described step (4), its concrete grammar is:

(4-1) current time is walked the location quantization function value LF of n ⁿbe mapped as the distribution density p of particle ⁿ(x);

(4-2) according to p ⁿx () calculates the population needing to increase/delete

(4-3) increase in grid cell x or delete individual particle, obtains the particle collection after upgrading with

(4-4) upgrade the speed of middle particle and position, then utilize them to detect and correct error, call sign distance field function phi ⁿ⁺¹;

(4-5) upgrade the speed of middle particle and position, obtain the spray particle collection of next time step n+1

In described step (5), concrete grammar comprises:

(5-1) marching cubes algorithm is used, from the symbolic distance field function φ representing flow surface ⁿ⁺¹middle extraction zero contour surface, obtains the triangular mesh surface of flow surface;

(5-2) structure particle local horizontal set algorithm is used, from particle collection the middle symbolic distance field function obtaining expression spray particle and then use marching cubes algorithm, obtain zero contour surface representing spray.

Beneficial effect of the present invention is:

(1) method based on geometric properties is adopted to realize the high precision modeling of convection cell animation, by obtaining geometric parameter, set up geometric model, geometric model and the process such as physical model is coupled, abundant excavation also uses the geometric properties such as the otherness of correlativity between fluid details and flow surface zones of different, and the robotization realizing high precision fluid animation generates;

(2) utilize location quantization function can effectively reflect the position that the details such as surface motions is rapid, spray of water produces and details number;

(3) the method has that precision is high, the feature of high efficiency.

Accompanying drawing explanation

Fig. 1 is schematic flow sheet of the present invention.

Embodiment:

Below in conjunction with accompanying drawing and embodiment, the invention will be further described.

As shown in Figure 1, the high precision fluid animation modeling method based on geometric properties mainly comprises following process:

Process 1: initialization;

The first step: initialization liquid surface.Make Ω represent computational fields, computational fields grid is designated as G.Make symbolic distance field function φ (x, t) express liquid surface, its functional value is the symbolic distance of t grid element center x to liquid surface.Obviously, on liquid surface, meet φ (x, t)=0, be designated as liquid internal meets φ (x, t) < 0, is designated as Ω ^-; Meet φ (x, t) > 0 in liquid external, be designated as Ω ⁺.

Second step: initialization fluid velocity field u.Specify the x, y, z component of three three-dimensional numerical values difference representation speed field u on computational fields grid G, at liquid external Ω ⁺on region, the value of three components is 0; And in liquid surface and inside nonzero value is appointed as according to concrete example in region.

3rd step: initialization detects and corrects particle collection P _c.In flow surface on inner and outside three-layer network lattice, initialization detects and corrects particle collection P _c.For P _cin each particle: its position is stochastic generation in inside and outside three-layer network lattice, and its radius is the size of 0.1-0.5 grid cell, and its speed obtains according to particle position interpolation from initialized velocity field u.

4th step: initialization represents the particle collection P of spray _sfor sky.

5th step: initialization boundary condition.The environmental model at place of being moved by fluid utilizes voxelization method to be converted into discernible border when equation calculates.It comprises two parts: solid and liquid boundary, air and liquid boundary.Use flag marker bit mark herein.Flag=0 represents the unit grid occupied by solid; Flag=1 express liquid unit; Flag=2 represents air element.

Process 2: according to current fluid velocity field u and particle collection P _sobtain the midrange speed field u in computational fields grid ^*.

The first step: make current time walk as n, the fluid velocity field in computational fields grid G is u ⁿ, SPH particle collection P _sin particle rapidity be in order to calculate dynamic fluid velocity field u ⁿ⁺¹, need the NS equation solving fluid motion

$\frac{1}{\mathrm{\&rho;}}\frac{\mathrm{D\&rho;}}{\mathrm{Dt}}+\&dtri;\&CenterDot;u=0---\left(1\right)$

$\frac{\&PartialD;u}{\&PartialD;t}=-u\&CenterDot;(\&dtri;u)-\frac{1}{\mathrm{\&rho;}}\&dtri;p+f---\left(2\right)$

Wherein, ρ is density, and u is fluid velocity, and p is pressure, and f is external force, and ▽ is divergence, and ▽ is gradient.

Second step: for each grid cell on region, its middle velocity field u ^*midrange speed field u is obtained by performing convection current and adding external force two steps ^*.Wherein, convection current step can use Semi-Lagrangian method to solve, namely add external force step and can be expressed as u ^*=u (x)+Δ tf (x, t), Δ t are time step.

3rd step: for Ω ⁺region need calculate P _su in each grid cell at particle collection place ^*.For this reason, external force is first added namely and then projected in grid G, namely wherein, ω _px () is for particle is to the contribution of grid

${\mathrm{\&omega;}}_p\left(x\right)=\left\{\begin{array}{cc}c(1-{\left|\right|x-x_p\left|\right|}^2/r_p^2)& \mathrm{if}{\left|\right|x-x_p\left|\right|}^2\&le;r_p^2\\ 0& \mathrm{else}\end{array}\right.$

Here, c is control constant, x _pfor particle position, r _pfor particle radii.

Process 3: according to midrange speed field, by solving the pressure Poisson equation of coupling Euler and Lagrange two kinds of models, obtains the fluid velocity field u in computational fields ⁿ⁺¹.

The first step: transformation equation (1) (2), obtain pressure Poisson equation.U is obtained according to process 2 method ^*after, equation (2) right side is remaining-1/ ρ ▽ p only, utilizes finite difference to launch,

u ⁿ⁺¹＝u ^*-▽pΔt/ρ (3)

Divergence is got to equation (3) both sides simultaneously, pressure Poisson equation can be obtained

▽·(▽pΔt/ρ)＝▽·u ^*-▽·u ⁿ⁺¹ (4)

Second step: the right scale of computing formula (4).In formula (4), ▽ u ⁿ⁺¹different in the diverse location value of computational fields.Wherein, region belongs to liquid internal, has strict incompressibility, therefore ▽ u ⁿ⁺¹=0.And at Ω ⁺region, liquid particles mixes with air, has certain compressibility, ▽ u ⁿ⁺¹≠ 0, its value can according to P in grid cell _sparticle number and intended particle number are estimated.According to mass-conservation equation (1), ρ can be expressed as ρ=N _pm _p/ V _p, wherein, Ν _pfor the P in grid cell _snumber of particles, m _pfor the quality of particle, V _pvolume shared by particle.Therefore equation (1) is variable is changed to

$\&dtri;\&CenterDot;u=-\frac{1}{N_p}\frac{{\mathrm{DN}}_p}{\mathrm{Dt}}---\left(5\right)$

In order to calculation stability, the segmentation integrated form herein in step service time calculate.Wherein Δ τ is the time step after segmentation, for the SPH particle number of n moment each grid cell, for the number of intended particle.

3rd step: solve pressure p.According to the first step and second step, the right side of equation (4) is known, then equation (4) is converted into the system of linear equations about pressure p.Use this system of linear equations of Conjugate Gradient Method With Preconditioning Efficient Solution can obtain pressure solution p.

4th step: obtain the fluid velocity field u in the computational fields G in n+1 moment ⁿ⁺¹.P is substituted into formula (3), velocity field new in fluid mass can be obtained.Finally by speed Extrapolation method, the speed in fluid mass is extrapolated in air-grid the velocity field u obtained in computational fields G ⁿ⁺¹.

Process 4: according to the velocity field u in obtained computational fields G ⁿ⁺¹, calculate the symbolic distance field function φ after upgrading ⁿ⁺¹.

The first step: the symbolic distance field function making current time walk n corresponding is φ ⁿ, utilize Semi-Lagrangian method to solve LS equation

$\frac{\&PartialD;\mathrm{\&phi;}}{\&PartialD;t}+u\&CenterDot;\&dtri;\mathrm{\&phi;}=0---\left(6\right)$

Obtain uncorrected symbolic distance field function here, u=u ⁿ⁺¹for the fluid velocity field after renewal.

Second step: use and detect correction particle collection P _cin particle pair carry out error-detecting and correction, call sign distance field function phi ⁿ⁺¹.Concrete, make s _pfor particle symbol, for the particle position of n+1 time step, if existed

$s_p{\mathrm{\&phi;}}_m^{n+1}\left(x_c^{n+1}\right)<0$

Then illustrate and produce particle escapes, namely detect that flow surface calculating exists error.Now, need to utilize escaped particles to define local LS function convection cell surface corrects.Timing, needs first to calculate the positive and negative particle collection E that escapes ⁺and E ^-the maximal value of the local LS of net point and minimum value then in the net point x of escaped particles collection place, two distance fields are merged.If | φ ⁺(x) |≤| φ ^-(x) |, then φ _c(x)=φ ⁺(x), otherwise φ _c(x)=φ ^-(x).By gained φ _cx () value replaces the distance field functional value of x position just the flow surface φ after error correction can be obtained ⁿ⁺¹.

Process 5: definition and calculating local axis LMA.

The first step: build envelope closed manifold, formed by sealing fluid interface or nonocclusive fluid interface and computational fields border.Envelope closed manifold border and around 3 layers of grid are called narrowband region.Arrowband is defined as wherein, for three-layer network lattice region outside liquid surface, for three-layer network lattice region inside liquid surface, express liquid surface.

Second step: definition local axis LMA is that above-mentioned closed stream morpheme is in arrowband Ω _bon the track in the inscribed circle center of circle.LMA is a point set.

3rd step: initialization LMA point set is empty.

3rd step: estimate the some p that LMA concentrates ^m.For this reason, the divergence value of the gradient of compute sign distance field function phi.Be expressed as

$D=\&dtri;\&CenterDot;(\&dtri;\mathrm{\&phi;})=\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;x}^2}+\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;y}^2}+\frac{{\&PartialD;}^2\mathrm{\&phi;}}{{\&PartialD;z}^2}=\frac{1}{\mathrm{\&Delta;}{h}^2}({\mathrm{\&phi;}}_{i-1,j,k}+{\mathrm{\&phi;}}_{i,j-1,k}+{\mathrm{\&phi;}}_{i,j,k-1}-6{\mathrm{\&phi;}}_{i,j,k}+{\mathrm{\&phi;}}_{i+1,j,k}+{\mathrm{\&phi;}}_{i,j+1,k}+{\mathrm{\&phi;}}_{i,j,k+1})$

Wherein, Δ h is spatial mesh size, is commonly defined as the size of space lattice.If namely the divergence value of grid cell is greater than certain threshold value, illustrates that this unit includes LMA point.Then the some p on LMA is similar to this grid element center point (net point) x ^m, be incorporated into existing LMA point and concentrate, be i.e. LMA=LMA ∪ x.

Process 6: definition location quantization function (Location and Quantization Function, LF).

The first step: definition interfaces on the distance of some x to LMA.Be expressed as:

$L(x,\mathrm{LMA})=\min (\underset{p^M\&Element;\mathrm{LMA}}{\inf }\left(\right||x-p^M|\left|\right),\mathrm{\&lambda;})---\left(7\right)$

Wherein, inf is floor value, and λ=10 Δ h is threshold value.

Second step: definition LF function.Make x be any point on computational fields, y is interface on a bit, then function LF is defined as:

$\mathrm{LF}\left(x\right)=\left\{\begin{array}{cc}L(x,\mathrm{LMA})& x\&Element;\&PartialD;\mathrm{\&Omega;}\\ k\left|\right|x-y\left|\right|+\mathrm{LF}\left(y\right)& x\&Element;{\mathrm{\&Omega;}}_b-\&PartialD;\mathrm{\&Omega;}\\ \&infin;& x\&Element;\mathrm{\&Omega;}-{\mathrm{\&Omega;}}_b\end{array}\right.---\left(8\right)$

Wherein, ∞ represents large constant, and k is controling parameters, the severe degree that in control inerface, the value of LF changes with distance.Y is interface the closest approach of upper distance x, L (x, LMA) represents the distance of some x to the LMA on interface.

Process 7: the LF functional value of net point on gauging surface.

The first step: build KD tree and store LMA point set.The KD tree of LMA build use the concentrated dimension with maximum span value of LMA point as cutting dimension, by the average in this dimension as cutting value.

Second step: the method using KD retrieval, the LF functional value of Fluid Computation net point on the surface.Can concentrate from LMA point fast from KD tree and obtain the nearest some p of distance x ^m.Finally, according to formula (7) gauging surface the LF value of upper net point x.

Process 8: calculate lF functional value on region.

The first step: initialization net point set: set up three set Far, Close and Accepted; Will the net point at place is inserted in Accepted, and being positioned at these net points on neighbours' net point insert in Close; Will on remaining net point insert in Far.

Second step: initialization LF value: in Accepted, the LF value of all net points is the input of algorithm, namely on LF value; During Close set and Far gather, the LF value of all net points is initialized as a large constant, i.e. LF=∞.

3rd step: set up with the LF value most rickle that is index, for preserving Close set.

4th step: the LF value upgrading all net point x in Close set: make Y={y} be the neighbours net point set of x in Accepted set, then the LF value of x is updated to

LF(x)＝min(LF(x),k||x-y||+LF(y))

5th step: upgrade the most rickle preserving Close collection according to the LF value that recalculates.

6th step: select most rickle top element T rial (namely there is in Close set the net point of minimum LF value), inserted Accepted and gather; Neighbours' point of all Far of the being arranged in set in Trail is transferred to Close set.

7th step: forward the 4th step to, iteration is run until most rickle be empty, turn the 8th step.

8th step: all in above-mentioned Accepted set on the LF value of net point be required by.

Process 9: calculate lF functional value on region.

The first step: initialization net point set: set up three set Far, Close and Accepted; Will the net point at place is inserted in Accepted, and being positioned at these net points on neighbours' net point insert in Close; Will on remaining net point insert in Far.

Second step: initialization LF value: in Accepted, the LF value of all net points is the input of algorithm, namely on LF value; During Close set and Far gather, the LF value of all net points is initialized as a large constant, i.e. LF=∞.

3rd step: set up with the LF value most rickle that is index, for preserving Close set.

4th step: the LF value upgrading all net point x in Close set: make Y={y} be the neighbours net point set of x in Accepted set, then the LF value of x is updated to

LF(x)＝min(LF(x),k||x-y||+LF(y))

5th step: upgrade the most rickle preserving Close collection according to the LF value that recalculates.

6th step: select most rickle top element T rial (namely there is in Close set the net point of minimum LF value), inserted Accepted and gather; Neighbours' point of all Far of the being arranged in set in Trail is transferred to Close set.

7th step: forward the 4th step to, iteration is run until most rickle be empty, turn the 8th step.

8th step: all in above-mentioned Accepted set on the LF value of net point be required by.

Process 10: calculate the LF functional value beyond narrow-band-domain.According to formula (8), be set to a large constant.

Process 11: the deformation extent value S on Fluid Computation surface.

The first step: definition list facial disfigurement function S.Flow surface warping function S is defined as the rate of change of the fluid volume that a period of time flow surface is surrounded

$S(t_i,t_{i-1})=\underset{t=t_{i-1},\mathrm{\&Delta;t}}{\overset{t_i}{\mathrm{\&Sigma;}}}\left(\right|A(t+\mathrm{\&Delta;t})-A\left(t\right)\left|\right)/A\left(t_{i-1}\right)---\left(9\right)$

Wherein, t _i-1and t _iresampling moment last time and current time respectively, Δ t is time step, and the fluid volume that A (t) is surrounded by fluid interface for t is expressed as

$A\left(t\right)=\underset{{\mathrm{\&Omega;}}_b}{\&Integral;}\mathrm{\&delta;}\left({\mathrm{\&phi;}}^t\right)|\&dtri;{\mathrm{\&phi;}}^t|\mathrm{dx}---\left(10\right)$

Wherein, δ function is Smeared-out Heaviside function derivative

$\mathrm{\&delta;}\left(\mathrm{\&phi;}\right)=\left\{\begin{array}{cc}0& \mathrm{\&phi;}<-\mathrm{\&epsiv;}\\ \frac{1}{2\mathrm{\&epsiv;}}+\frac{1}{2\mathrm{\&epsiv;}}\cos \left(\frac{\mathrm{\&pi;\&phi;}}{\mathrm{\&epsiv;}}\right)& -\mathrm{\&epsiv;}\&le;\mathrm{\&phi;}\&le;\mathrm{\&epsiv;}\\ 0& \mathrm{\&epsiv;}<\mathrm{\&phi;}\end{array}\right..$

Here, ε=1.5 Δ h, Δ h are spatial mesh size.

Second step: the S value of calculating.In the calculation, formula (10) can be converted into

$A\left(t\right)=\underset{(i,j,k)\&Element;{\mathrm{\&Omega;}}_b}{\mathrm{\&Sigma;}}\frac{1}{\mathrm{\&Delta;h}}\mathrm{\&delta;}\left(\mathrm{\&phi;}(i,j,k,t)\right)\sqrt{\mathrm{\&Phi;}{\left(I\right)}^2+\mathrm{\&Phi;}{\left(J\right)}^2+\mathrm{\&Phi;}{\left(K\right)}^2}$

Φ(I)＝φ(i+1,j,k,t)-φ(i,j,k,t)，

Φ(J)＝φ(i,j+1,k,t)-φ(i,j,k,t)，

Φ(K)＝φ(i,j,k+1,t)-φ(i,j,k,t)，

Wherein (i, j, k) is arrowband Ω _bon net point.Like this, according to process 4 symbolic distance field function φ ⁿ⁺¹the parameter obtained, brings formula (9) (10) into and can calculate surface deformation degree S.

Process 12: upgrade LF functional value.

The first step: determine the character that current time walks.If the value of S is less than threshold value θ, instruction card facial disfigurement is little, and without the need to recalculating geometric properties, this time step belongs to convection current step; Otherwise this time step belongs to and reinitializes step.

Second step: if convection current step, then only as flow surface, convection current need be carried out to it and upgrade.The now renewal of LF is expressed as

$\frac{\&PartialD;F}{\&PartialD;t}=-u\&CenterDot;\&dtri;F---\left(11\right)$

Wherein, F=LF (x, t), u=u ⁿ⁺¹for the fluid velocity field in the grid G that current time has obtained.Like this according to the LF value LF of current time ⁿwith the G medium velocity field u after renewal ⁿ⁺¹, utilize Semi-Lagrangian method solution formula (11), the LF value LF after renewal can be obtained ⁿ⁺¹.

3rd step: if reinitialize step, then need to recalculate LF functional value according to the method for process (7)-process (10).

Process 13: according to the fundamental function obtained from geometric model calculating, increase or delete the particle of particle set.

The first step: LF functional value LF current time being walked n ⁿbe mapped as the distribution density p of particle ⁿ(x).Because LF is less at the functional value of deep camber, elongated area and borderline region, and particle needed for these regions is more on the contrary, therefore, first converts LF

$p^n\left(x\right)=1.0-\frac{{\mathrm{LF}}^n\left(x\right)-\min \left({\mathrm{LF}}^n\right)}{\max \left({\mathrm{LF}}^n\right)-\min \left({\mathrm{LF}}^n\right)}.$

Like this, grid (i.e. deep camber, elongated area and the interfacial boundary etc.) p that LF value is less ⁿx () value is comparatively large, thus meet the demand in these areal distribution more multiparticle.Then by normalization, the distribution density of particle is obtained

$p^n\left(x\right)=\frac{p^n\left(x\right)}{\underset{x}{\mathrm{\&Sigma;}}{p}^n\left(x\right)}$

Second step: according to p ⁿx () calculates the population needing to increase/delete for

$N_{\mathrm{del}}^n\left(x\right)=\mathrm{\&kappa;}(\underset{x}{\mathrm{\&Sigma;}}{N}_p^n\left(x\right)\&times;p^n\left(x\right)-N_p^n\left(x\right))$

Wherein, for the population in current grid unit x, for the total number of particle in computational fields, κ is controling parameters, controls the severe degree of particle number change.

3rd step: increase in grid cell x or delete individual particle, obtains the particle collection after upgrading with in each grid cell x, if then random erasure in grid individual particle.Otherwise, stochastic generation individual particle, particle position produces in grid cell x internal random, and speed is initialized as zero.

4th step: determine the classification of the new particle under current time step n and insert respectively with in insert with in with obtain with for new particle p, if namely this particle is positioned at region, then belong to its symbol is negative.Otherwise then this particle is positioned at region.Here, for the uncorrected symbolic distance field function obtained in process (4).If particle is positioned at interface ω layer (ω=3) grid in, namely then this particle belongs to its symbol is just.If then particle belongs to

Process 14: upgrade the speed of middle particle and position, right implement error-detecting and correct acquisition symbolic distance field function φ ⁿ⁺¹.

The first step: upgrade the speed of middle particle.Right in particle, its speed the method of trilinear interpolation can be utilized from G medium velocity field u ⁿ⁺¹obtain.Namely wherein for the position of particle.

Second step: upgrade the position of middle particle.With as input, pass through solving equation the particle position after renewal can be obtained if in the air-grid of negative particle escapes outside liquid surface three-layer network lattice in, then by its from in remove, add to in.

3rd step: obtain particle collection new in n+1 time step above-mentioned in after all particle rapidities and location updating, particle collection newly

4th step: utilize new particle collection detect and correct obtain due to for the uncorrected symbolic distance field function obtained in process (4).Due to particle collection in the present invention distribution according to geometric properties calculate gained, therefore it is right in error can detect more accurately and correct more accurately, improve the precision of flow surface.Concrete, make s _pfor particle symbol, if existed

$s_p{\mathrm{\&phi;}}_m^{n+1}\left(x_c^{n+1}\right)<0;$

Then particle escapes is described, namely detects that flow surface calculating exists error.Now, need to utilize escaped particles to define local LS function convection cell surface corrects.Timing, needs first to calculate the positive and negative particle collection E that escapes ⁺and E ^-the maximal value of the local LS of net point and minimum value then in the net point x of escaped particles collection place, two distance fields are merged.If | φ ⁺(x) |≤| φ ^-(x) |, then φ _c(x)=φ ⁺(x), otherwise φ _c(x)=φ ^-(x).By gained φ _cx () value replaces the distance field functional value of x position just the flow surface φ after error correction can be obtained ⁿ⁺¹.

Process 15: upgrade the speed of middle particle and position, obtain the spray particle collection of next time step n+1

The first step: upgrade the speed of middle particle due to u ⁿ⁺¹calculating be particle rapidity is projected in grid and solve grid and particle coupling pressure equation obtain.Therefore, only need by u ⁿ⁺¹back projection to particle can obtain new particle rapidity.Concrete, first midrange speed field u ^*to u ⁿ⁺¹increment back projection on particle, i.e. Δ u _p=Σ (u ⁿ⁺¹-u ^*) ω _px (), then the particle rapidity after upgrading is parameter s represents that particle rapidity is subject to the influence degree of pressure.

Second step: upgrade middle particle position according to solving equation particle position after can upgrading if within Particles Moving to flow surface, namely then from in remove this particle.

3rd step: obtain particle collection new in n+1 time step above-mentioned after middle particle rapidity and location updating, obtain particle collection newly

Process 16: from above-mentioned result of calculation φ ⁿ⁺¹with middle triangular mesh surface of extracting expression flow surface and fluid spray, and can fluid animation be generated after Realistic Rendering is carried out to it.

The first step: use marching cubes algorithm, from the symbolic distance field function φ representing flow surface ⁿ⁺¹middle extraction zero contour surface, obtains the triangular mesh surface of flow surface.

Second step: use and build particle local horizontal set algorithm, from particle collection the middle symbolic distance field function obtaining expression spray particle be defined as:

${\mathrm{\&phi;}}_p^{n+1}\left(x\right)=|x-\stackrel{\&OverBar;}{x}|-\stackrel{\&OverBar;}{r}$

Wherein, for the mean place of particle periphery neighbor particle, be designated as for the mean radius of particle periphery neighbor particle, be designated as here i is the neighbor particle label of current particle.W _ifor weighted value.

$w_i=\frac{k\left(\right|x-x_i|/R)}{\underset{i}{\mathrm{\&Sigma;}}k\left(\right|x-x_j|/R)}$

Here, k is that kernel function is defined as k (d)=max (0, (1-d ²) ³).

After obtaining the symbolic distance field function of spray particle, re-use marching cubes algorithm, obtain zero contour surface, this contour surface is the triangular mesh surface of spray.

3rd step: use global illumination algorithm, carries out Realistic Rendering to the triangular mesh surface of the flow surface obtained and fluid spray, can obtain fluid animation.

By reference to the accompanying drawings the specific embodiment of the present invention is described although above-mentioned; but not limiting the scope of the invention; one of ordinary skill in the art should be understood that; on the basis of technical scheme of the present invention, those skilled in the art do not need to pay various amendment or distortion that creative work can make still within protection scope of the present invention.

Claims (7)

1., based on a high precision fluid animation modeling method for geometric properties, it is characterized in that: comprise the following steps:

(1) starting condition in convection cell physical model and boundary condition carry out initialization;

(2) from the physical model of fluid, obtain the fluid velocity field in the midrange speed field required for geometric properties calculating, current time step and symbolic distance field function parameter;

(3) build the geometric model of fluid animation, robotization location, quantification and tracking are carried out in the details sensitizing range in convection cell animation;

(4) according to the fundamental function of geometric properties, the particle collection of physical model is controlled, coupling geometric model and physical model;

(5) extract the triangular mesh surface representing flow surface and fluid spray, and Realistic Rendering is carried out to it, generate a frame fluid animation;

(6) check the making whether completing whole fluid animation, if do not completed, then return step (2), proceed to make.

2. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 1, it is characterized in that: in described step (1), the initialized step of physical model comprises:

(1-1) the fluid velocity field u in initial background grid;

(1-2) initialization represents the symbolic distance field function φ of flow surface;

(1-3) initialization convection cell surface carries out detecting the particle collection P corrected _cwith the particle collection P representing fluid spray _s;

(1-4) initialization border, namely uses marker bit to represent the type of each grid cell in fluid calculation territory; Need, by solid obstacle by computational fields grid cell size discrete, to be set to 1 by the grid cell marker bit that solid occupies after discretize, otherwise to be set to 0 for this reason.

3. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 1, is characterized in that: in described step (2), the concrete grammar obtaining geometric properties from the physical model of fluid is as follows:

(2-1) according to the fluid velocity field u in initial background grid, particle collection P _sin speed calculate midrange speed field in background grid;

(2-2) according to midrange speed field, by solving the pressure Poisson equation of coupling Euler and Lagrange two kinds of models, the fluid velocity field u in current time step is obtained ⁿ⁺¹;

(2-3) according to obtained current time step fluid velocity field, the symbolic distance field function φ after upgrading is calculated ⁿ⁺¹.

4. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 1, it is characterized in that: in described step (3), set up the geometric model of fluid animation, the details sensitizing range in convection cell animation carry out robotization location, quantification and tracking basic step as follows:

(3-1) build envelope closed manifold, is formed by sealing fluid interface or nonocclusive fluid interface and computational fields border, seal closed manifold border and surrounding 3 layers of grid are called narrowband region;

(3-2) definition local axis is for being positioned at arrowband Ω in envelope closed manifold _bon the track in the inscribed circle center of circle;

(3-3) according to the symbolic distance field function φ obtained in Step2 ⁿ⁺¹, use divergence value determining method to calculate LMA;

(3-4) according to local axis, definition location quantization function;

(3-5) according to the definition of location quantization function, by using the method for KD retrieval, the location quantization function value of Fluid Computation net point on the surface;

(3-6) according to the definition of location quantization function, by using the location quantization function value extrapolation of fast marching algorithms convection current shape frontier point, the location quantization function value of other net points in narrow-band-domain is calculated;

(3-7) according to the definition of location quantization function, the location quantization function value in setting background grid beyond narrow-band-domain is large constant, is appointed as N _max× Δ h, wherein N _maxfor computational fields is in the maximal value of x, y, z three axially meshes number, Δ h is the size of grid cell;

(3-8) deformation extent on accumulation Fluid Computation surface in each time step, if deformation extent is less than setting threshold value, then needs to carry out convection current renewal to location quantization function; Otherwise need to reinitialize location quantization function.

5. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 4, it is characterized in that: in described step (3-4), concrete grammar is: according to local axis, definition location quantization function, for stream shape borderline net point, its functional value is given threshold value λ and the smaller value of this point to locally half-breadth value; In narrowband region, the location quantization function value of other net points is that this point adds the location quantization function value of shape border, upper reaches corresponding point to the bee-line flowing shape border; Location quantization function value for other net points on computational fields is set to a large constant, and location quantization function can realize smooth location to deep camber, elongated and edge details sensitizing range and quantification.

6. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 1, is characterized in that: in described step (4), its concrete grammar is:

(4-1) current time is walked the location quantization function value LF of n ⁿbe mapped as the distribution density p of particle ⁿ(x);

(4-2) according to p ⁿx () calculates the population needing to increase/delete

(4-3) increase in grid cell x or delete individual particle, obtains the particle collection after upgrading with

(4-4) upgrade the speed of middle particle and position, then utilize them to detect and correct error, call sign distance field function phi ⁿ⁺¹;

(4-5) upgrade the speed of middle particle and position, obtain the spray particle collection of next time step n+1

7. a kind of high precision fluid animation modeling method based on geometric properties as claimed in claim 1, is characterized in that: in described step (5), concrete grammar comprises:

(5-1) marching cubes algorithm is used, from the symbolic distance field function φ representing flow surface ⁿ⁺¹middle extraction zero contour surface, obtains the triangular mesh surface of flow surface;

(5-2) structure particle local horizontal set algorithm is used, from particle collection the middle symbolic distance field function obtaining expression spray particle and then use marching cubes algorithm, obtain zero contour surface representing spray.