Method matrixFreeRHS

SPlisHSPlasH/PF/TimeStepPF.cpp:371–593

\brief compute the right hand side of the system in a matrix-free fashion and store the result in result*/

/** \brief compute the right hand side of the system in a matrix-free fashion and store the result in result*/
void SPH::TimeStepPF::matrixFreeRHS(const VectorXr & x, VectorXr & result)
{
	Simulation *sim = Simulation::getCurrent();
	const Real h = TimeManager::getCurrent()->getTimeStepSize();;
	const unsigned int nFluids = sim->numberOfFluidModels();
	const unsigned int nBoundaries = sim->numberOfBoundaryModels();

	//////////////////////////////////////////////////////////////////////////
	// helper functions
	//////////////////////////////////////////////////////////////////////////
	// constraint value
	const auto calculateC = [&](const unsigned int fluidModelIndex, const unsigned int i, std::vector<Vector3r> & p) -> Real
	{
		const FluidModel * model = sim->getFluidModel(fluidModelIndex);
		// Compute current density for particle i
		Real density = model->getVolume(i) * sim->W_zero();
		const Vector3r &xi = p[0];
		unsigned int counter = 1;

		for (unsigned int pid = 0; pid < nFluids; pid++)
		{
			const FluidModel *fm_neighbor = sim->getFluidModelFromPointSet(pid);
			for (unsigned int j = 0; j < sim->numberOfNeighbors(fluidModelIndex, pid, i); j++)
			{
				const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
				const Vector3r & xj = p[counter++];
				density += fm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
			}
		}
		// influence of boundary on density
		if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)
		{
			forall_boundary_neighbors(
				// Boundary: Akinci2012
				density += bm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
			);
		}
		else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017)
		{
			forall_density_maps(
				density += rho;
			);
		}
		else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019)
		{
			forall_volume_maps(
				density += Vj * sim->W(xi - xj);
			);
		}
		// constraint value = density / density0 - 1
		const auto C = density - 1;
		// pressure clamping
		return (C < 0) ? 0 : C;
	};
	// constraint gradient
	const auto calculateNablaC = [&](const unsigned int fluidModelIndex, const unsigned int i, std::vector<Vector3r> & p) -> std::vector<Vector3r>
	{
		std::vector<Vector3r> nablaC(p.size());