Method step

Tools/VolumeSampling/SPHVolumeSampling.cpp:112–211

void SPHVolumeSampling::step(Real& avg_density_err)
{
	avg_density_err = 0.0;

	//////////////////////////////////////////////////////////////////////////
	// SPH optimization of sampling points
	//////////////////////////////////////////////////////////////////////////
	
	const int numParticles = (int)m_x.size();	
	Real V0 = m_diameter * m_diameter * m_diameter;

	const Real diameter2 = m_diameter* m_diameter;
	const Real supportRadius = m_radius* static_cast<Real>(4.0);

	START_TIMING("neighborhoodSearch");
	m_neighborhoodSearch->find_neighbors();
	STOP_TIMING_AVG

	computeDensities(m_densities, m_mass);

	// Compute cohesion => Becker 2007
	#pragma omp parallel default(shared)
	{
		#pragma omp for schedule(static)  
		for (int i = 0; i < numParticles; i++)
		{
			const Vector3r &xi = m_x[i];
			Vector3r &corr = m_corrs[i];
			corr.setZero();

 			for (unsigned int j = 0; j < numberOfNeighbors(0, i); j++)
 			{
 				const unsigned int neighborIndex = getNeighbor(0, i, j);
 				const Vector3r &xj = m_x[neighborIndex];
 
 				const Vector3r xixj = xi - xj;
 				const Real r2 = xixj.dot(xixj);
				if (r2 > 1e-6)
				{
					const Real C = (xi - xj).norm() - m_diameter;
					Vector3r xixj_n = xixj;
					xixj_n.normalize();
					corr -= m_cohesion * m_mass/ m_densities[neighborIndex] * C * xixj_n * m_kernelFct(xi - xj);
				}
 			}

			if (m_adhesion != 0.0)
			{
				// Collision detection
				Vector3r cp, normal;
				const double dist = distance(xi, 0.0, normal, cp);
				if (dist < supportRadius)
				{
					corr -= m_adhesion * m_mass / m_densities[i] * (xi - cp) * m_kernelFct(xi - cp);
				}
			}
		}
	}