1+ #include " FluidModel.h"
2+ #include " PositionBasedDynamics/PositionBasedDynamics.h"
3+ #include " PositionBasedDynamics/SPHKernels.h"
4+
5+ using namespace PBD ;
6+
7+ FluidModel::FluidModel () :
8+ m_particles()
9+ {
10+ m_density0 = 1000 .0f ;
11+ m_particleRadius = 0 .025f ;
12+ viscosity = 0 .02f ;
13+ m_neighborhoodSearch = NULL ;
14+ }
15+
16+ FluidModel::~FluidModel (void )
17+ {
18+ cleanupModel ();
19+ }
20+
21+ void FluidModel::cleanupModel ()
22+ {
23+ m_particles.release ();
24+ m_lambda.clear ();
25+ m_density.clear ();
26+ m_deltaX.clear ();
27+ delete m_neighborhoodSearch;
28+ }
29+
30+ void FluidModel::reset ()
31+ {
32+ const unsigned int nPoints = m_particles.size ();
33+
34+ for (unsigned int i=0 ; i < nPoints; i++)
35+ {
36+ const Eigen::Vector3f& x0 = m_particles.getPosition0 (i);
37+ m_particles.getPosition (i) = x0;
38+ m_particles.getLastPosition (i) = m_particles.getPosition (i);
39+ m_particles.getOldPosition (i) = m_particles.getPosition (i);
40+ m_particles.getVelocity (i).setZero ();
41+ m_particles.getAcceleration (i).setZero ();
42+ m_deltaX[i].setZero ();
43+ m_lambda[i] = 0 .0f ;
44+ m_density[i] = 0 .0f ;
45+ }
46+ }
47+
48+ ParticleData & PBD::FluidModel::getParticles ()
49+ {
50+ return m_particles;
51+ }
52+
53+ void FluidModel::initMasses ()
54+ {
55+ const int nParticles = (int ) m_particles.size ();
56+ const float diam = 2 .0f *m_particleRadius;
57+
58+ #pragma omp parallel default(shared)
59+ {
60+ #pragma omp for schedule(static)
61+ for (int i = 0 ; i < nParticles; i++)
62+ {
63+ m_particles.setMass (i, 0 .8f * diam*diam*diam * m_density0); // each particle represents a cube with a side length of r
64+ // mass is slightly reduced to prevent pressure at the beginning of the simulation
65+ }
66+ }
67+ }
68+
69+
70+ /* * Resize the arrays containing the particle data.
71+ */
72+ void FluidModel::resizeFluidParticles (const unsigned int newSize)
73+ {
74+ m_particles.resize (newSize);
75+ m_lambda.resize (newSize);
76+ m_density.resize (newSize);
77+ m_deltaX.resize (newSize);
78+ }
79+
80+
81+ /* * Release the arrays containing the particle data.
82+ */
83+ void FluidModel::releaseFluidParticles ()
84+ {
85+ m_particles.release ();
86+ m_lambda.clear ();
87+ m_density.clear ();
88+ m_deltaX.clear ();
89+ }
90+
91+ void FluidModel::initModel (const unsigned int nFluidParticles, Eigen::Vector3f* fluidParticles, const unsigned int nBoundaryParticles, Eigen::Vector3f* boundaryParticles)
92+ {
93+ releaseFluidParticles ();
94+ resizeFluidParticles (nFluidParticles);
95+
96+ // init kernel
97+ CubicKernel::setRadius (m_supportRadius);
98+
99+ // copy fluid positions
100+ #pragma omp parallel default(shared)
101+ {
102+ #pragma omp for schedule(static)
103+ for (int i = 0 ; i < (int )nFluidParticles; i++)
104+ {
105+ m_particles.getPosition0 (i) = fluidParticles[i];
106+ }
107+ }
108+
109+ m_boundaryX.resize (nBoundaryParticles);
110+ m_boundaryPsi.resize (nBoundaryParticles);
111+
112+ // copy boundary positions
113+ #pragma omp parallel default(shared)
114+ {
115+ #pragma omp for schedule(static)
116+ for (int i = 0 ; i < (int )nBoundaryParticles; i++)
117+ {
118+ m_boundaryX[i] = boundaryParticles[i];
119+ }
120+ }
121+
122+ // initialize masses
123+ initMasses ();
124+
125+ // ////////////////////////////////////////////////////////////////////////
126+ // Compute value psi for boundary particles (boundary handling)
127+ // (see Akinci et al. "Versatile rigid - fluid coupling for incompressible SPH", Siggraph 2012
128+ // ////////////////////////////////////////////////////////////////////////
129+
130+ // Search boundary neighborhood
131+ NeighborhoodSearchSpatialHashing neighborhoodSearchSH (nBoundaryParticles, m_supportRadius);
132+ neighborhoodSearchSH.neighborhoodSearch (&m_boundaryX[0 ]);
133+
134+ unsigned int **neighbors = neighborhoodSearchSH.getNeighbors ();
135+ unsigned int *numNeighbors = neighborhoodSearchSH.getNumNeighbors ();
136+
137+ #pragma omp parallel default(shared)
138+ {
139+ #pragma omp for schedule(static)
140+ for (int i = 0 ; i < (int ) nBoundaryParticles; i++)
141+ {
142+ float delta = CubicKernel::W_zero ();
143+ for (unsigned int j = 0 ; j < numNeighbors[i]; j++)
144+ {
145+ const unsigned int neighborIndex = neighbors[i][j];
146+ delta += CubicKernel::W (m_boundaryX[i] - m_boundaryX[neighborIndex]);
147+ }
148+ const float volume = 1 .0f / delta;
149+ m_boundaryPsi[i] = m_density0 * volume;
150+ }
151+ }
152+
153+
154+ // Initialize neighborhood search
155+ if (m_neighborhoodSearch == NULL )
156+ m_neighborhoodSearch = new NeighborhoodSearchSpatialHashing (m_particles.size (), m_supportRadius);
157+ m_neighborhoodSearch->setRadius (m_supportRadius);
158+
159+ reset ();
160+ }
0 commit comments