Strong Fluid-Rigid Coupling

SPlisHSPlasH/BoundaryModel_Akinci2012.cpp

@@ -6,6 +6,7 @@
 #include "Utilities/Logger.h"
 #include "NeighborhoodSearch.h"
 #include "Simulation.h"
+#include "StrongCouplingBoundarySolver.h"
 
 using namespace SPH;
 
@@ -18,6 +19,11 @@ BoundaryModel_Akinci2012::BoundaryModel_Akinci2012() :
 {		
 	m_sorted = false;
 	m_pointSetIndex = 0;
+
+	addField({ "position", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &getPosition(i)[0]; }, true });
+	addField({ "position0", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &getPosition0(i)[0]; } });
+	addField({ "velocity", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &getVelocity(i)[0]; }, true });
+	addField({ "volume", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &getVolume(i); }, true });
 }
 
 BoundaryModel_Akinci2012::~BoundaryModel_Akinci2012(void)
@@ -36,7 +42,6 @@ void BoundaryModel_Akinci2012::reset()
 	// positions and velocities are already updated by updateBoundaryParticles
 	if (!m_rigidBody->isDynamic() && !m_rigidBody->isAnimated())
 	{
-		// reset velocities and accelerations
 		for (int j = 0; j < (int)numberOfParticles(); j++)
 		{
 			m_x[j] = m_x0[j];
@@ -45,27 +50,61 @@ void BoundaryModel_Akinci2012::reset()
 	}
 }
 
+void SPH::BoundaryModel_Akinci2012::addField(const FieldDescription& field) {
+	m_fields.push_back(field);
+	std::sort(m_fields.begin(), m_fields.end(), [](FieldDescription& i, FieldDescription& j) -> bool { return (i.name < j.name); });
+}
+
+const FieldDescription& SPH::BoundaryModel_Akinci2012::getField(const std::string& name) {
+	unsigned int index = 0;
+	for (auto i = 0; i < m_fields.size(); i++) {
+		if (m_fields[i].name == name) {
+			index = i;
+			break;
+		}
+	}
+	return m_fields[index];
+}
+
+void SPH::BoundaryModel_Akinci2012::removeFieldByName(const std::string& fieldName) {
+	for (auto it = m_fields.begin(); it != m_fields.end(); it++) {
+		if (it->name == fieldName) {
+			m_fields.erase(it);
+			break;
+		}
+	}
+}
+
 void BoundaryModel_Akinci2012::computeBoundaryVolume()
 {
 	Simulation *sim = Simulation::getCurrent();
 	const unsigned int nFluids = sim->numberOfFluidModels();
 	NeighborhoodSearch *neighborhoodSearch = Simulation::getCurrent()->getNeighborhoodSearch();
-
+	StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent();
 	const unsigned int numBoundaryParticles = numberOfParticles();
 
 	#pragma omp parallel default(shared)
 	{
 		#pragma omp for schedule(static)  
 		for (int i = 0; i < (int)numBoundaryParticles; i++)
 		{
-			Real delta = sim->W_zero();
+			Real delta;
+			if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) {
+				delta = bs->W_zero();
+			} else {
+				delta = sim->W_zero();
+			}
 			for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++)
 			{
 				BoundaryModel_Akinci2012 *bm_neighbor = static_cast<BoundaryModel_Akinci2012*>(sim->getBoundaryModelFromPointSet(pid));
 				for (unsigned int j = 0; j < neighborhoodSearch->point_set(m_pointSetIndex).n_neighbors(pid, i); j++)
 				{
 					const unsigned int neighborIndex = neighborhoodSearch->point_set(m_pointSetIndex).neighbor(pid, i, j);
-					delta += sim->W(getPosition(i) - bm_neighbor->getPosition(neighborIndex));
+					if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) {
+						delta += bs->W(getPosition(i) - bm_neighbor->getPosition(neighborIndex));
+					} else {
+						delta += sim->W(getPosition(i) - bm_neighbor->getPosition(neighborIndex));
+					}
 				}
 			}
 			const Real volume = static_cast<Real>(1.0) / delta;

SPlisHSPlasH/BoundaryModel_Akinci2012.h

@@ -6,6 +6,7 @@
 
 #include "BoundaryModel.h"
 #include "SPHKernels.h"
+#include "FluidModel.h"
 
 
 namespace SPH 
@@ -33,8 +34,16 @@ namespace SPH
 			std::vector<Vector3r> m_x;
 			std::vector<Vector3r> m_v;
 			std::vector<Real> m_V;
+			std::vector<FieldDescription> m_fields;
 
 		public:
+			void addField(const FieldDescription& field);
+			const std::vector<FieldDescription>& getFields() {return m_fields;}
+			const FieldDescription& getField(const unsigned int i) {return m_fields[i];}
+			const FieldDescription& getField(const std::string& name);
+			const unsigned int numberOfFields() {return static_cast<unsigned int>(m_fields.size());}
+			void removeFieldByName(const std::string& fieldName);
+
 			unsigned int numberOfParticles() const { return static_cast<unsigned int>(m_x.size()); }
 			unsigned int getPointSetIndex() const { return m_pointSetIndex; }
 			bool isSorted() const { return m_sorted; }
@@ -50,7 +59,8 @@ namespace SPH
 			virtual void loadState(BinaryFileReader &binReader);
 
 			void initModel(RigidBodyObject *rbo, const unsigned int numBoundaryParticles, Vector3r *boundaryParticles);
-
+
+
 			FORCE_INLINE Vector3r &getPosition0(const unsigned int i)
 			{
 				return m_x0[i];

SPlisHSPlasH/CMakeLists.txt

@@ -209,6 +209,7 @@ add_library(SPlisHSPlasH
 	BoundaryModel_Bender2019.h
 	BoundaryModel_Koschier2017.cpp
 	BoundaryModel_Koschier2017.h
+	DynamicRigidBody.h
 	Emitter.cpp
 	Emitter.h
 	EmitterSystem.cpp
@@ -232,6 +233,8 @@ add_library(SPlisHSPlasH
 	NonPressureForceBase.h
 	XSPH.cpp
 	XSPH.h
+	StrongCouplingBoundarySolver.cpp
+	StrongCouplingBoundarySolver.h
 
 
 	${WCSPH_HEADER_FILES}

SPlisHSPlasH/DFSPH/SimulationDataDFSPH.cpp

@@ -9,7 +9,8 @@ SimulationDataDFSPH::SimulationDataDFSPH() :
 	m_pressure_rho2(),
 	m_pressure_rho2_V(),
 	m_pressureAccel(),
-	m_density_adv()
+	m_density_adv(),
+	m_density_adv_no_boundary()
 {
 }
 
@@ -29,6 +30,7 @@ void SimulationDataDFSPH::init()
 	m_pressure_rho2.resize(nModels);
 	m_pressure_rho2_V.resize(nModels);
 	m_pressureAccel.resize(nModels);
+	m_density_adv_no_boundary.resize(nModels);
 
 	for (unsigned int i = 0; i < nModels; i++)
 	{
@@ -38,6 +40,7 @@ void SimulationDataDFSPH::init()
 		m_pressure_rho2[i].resize(fm->numParticles(), 0.0);
 		m_pressure_rho2_V[i].resize(fm->numParticles(), 0.0);
 		m_pressureAccel[i].resize(fm->numParticles(), Vector3r::Zero());
+		m_density_adv_no_boundary[i].resize(fm->numParticles(), 0.0);
 	}
 }
 
@@ -53,12 +56,14 @@ void SimulationDataDFSPH::cleanup()
 		m_pressure_rho2[i].clear();
 		m_pressure_rho2_V[i].clear();
 		m_pressureAccel[i].clear();
+		m_density_adv_no_boundary[i].clear();
 	}
 	m_factor.clear();
 	m_density_adv.clear();
 	m_pressure_rho2.clear();
 	m_pressure_rho2_V.clear();
 	m_pressureAccel.clear();
+	m_density_adv_no_boundary.clear();
 }
 
 void SimulationDataDFSPH::reset()
@@ -76,6 +81,7 @@ void SimulationDataDFSPH::reset()
 			m_pressure_rho2_V[i][j] = 0.0;
 			m_factor[i][j] = 0.0;
 			m_pressureAccel[i][j].setZero();
+			m_density_adv_no_boundary[i][j] = 0.0;
 		}
 	}
 }

SPlisHSPlasH/DFSPH/SimulationDataDFSPH.h

@@ -25,6 +25,8 @@ namespace SPH
 			std::vector<std::vector<Real>> m_factor;
 			/** \brief advected density */
 			std::vector<std::vector<Real>> m_density_adv;
+			/** \brief for strong coupling (Gissler2019) */
+			std::vector<std::vector<Real>> m_density_adv_no_boundary;
 
 			/** \brief stores \f$\frac{p}{\rho^2}\f$ value of the constant density solver */
 			std::vector<std::vector<Real>> m_pressure_rho2;
@@ -85,6 +87,18 @@ namespace SPH
 				m_density_adv[fluidIndex][i] = d;
 			}
 
+			FORCE_INLINE const Real getDensityAdvNoBoundary(const unsigned int fluidIndex, const unsigned int i) const {
+				return m_density_adv_no_boundary[fluidIndex][i];
+			}
+
+			FORCE_INLINE Real& getDensityAdvNoBoundary(const unsigned int fluidIndex, const unsigned int i) {
+				return m_density_adv_no_boundary[fluidIndex][i];
+			}
+
+			FORCE_INLINE void setDensityAdvNoBoundary(const unsigned int fluidIndex, const unsigned int i, const Real d) {
+				m_density_adv_no_boundary[fluidIndex][i] = d;
+			}
+
 			FORCE_INLINE const Real getPressureRho2(const unsigned int fluidIndex, const unsigned int i) const
 			{
 				return m_pressure_rho2[fluidIndex][i];