diff --git a/SPlisHSPlasH/BoundaryModel_Akinci2012.cpp b/SPlisHSPlasH/BoundaryModel_Akinci2012.cpp index 46ec925b..ffd0a43f 100644 --- a/SPlisHSPlasH/BoundaryModel_Akinci2012.cpp +++ b/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,12 +50,37 @@ 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) @@ -58,14 +88,23 @@ void BoundaryModel_Akinci2012::computeBoundaryVolume() #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(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(1.0) / delta; diff --git a/SPlisHSPlasH/BoundaryModel_Akinci2012.h b/SPlisHSPlasH/BoundaryModel_Akinci2012.h index d0ee744b..a40e1bc6 100644 --- a/SPlisHSPlasH/BoundaryModel_Akinci2012.h +++ b/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 m_x; std::vector m_v; std::vector m_V; + std::vector m_fields; public: + void addField(const FieldDescription& field); + const std::vector& 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(m_fields.size());} + void removeFieldByName(const std::string& fieldName); + unsigned int numberOfParticles() const { return static_cast(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]; diff --git a/SPlisHSPlasH/CMakeLists.txt b/SPlisHSPlasH/CMakeLists.txt index b0238866..a38ec96e 100644 --- a/SPlisHSPlasH/CMakeLists.txt +++ b/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} diff --git a/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.cpp b/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.cpp index 25fe8855..f0fce750 100644 --- a/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.cpp +++ b/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; } } } diff --git a/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.h b/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.h index e99d8a5f..8fb4c880 100644 --- a/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.h +++ b/SPlisHSPlasH/DFSPH/SimulationDataDFSPH.h @@ -25,6 +25,8 @@ namespace SPH std::vector> m_factor; /** \brief advected density */ std::vector> m_density_adv; + /** \brief for strong coupling (Gissler2019) */ + std::vector> m_density_adv_no_boundary; /** \brief stores \f$\frac{p}{\rho^2}\f$ value of the constant density solver */ std::vector> 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]; diff --git a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp index afd9cb04..bf04ed2a 100644 --- a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp +++ b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp @@ -9,6 +9,9 @@ #include "SPlisHSPlasH/BoundaryModel_Akinci2012.h" #include "SPlisHSPlasH/BoundaryModel_Koschier2017.h" #include "SPlisHSPlasH/BoundaryModel_Bender2019.h" +#include "SPlisHSPlasH/StrongCouplingBoundarySolver.h" +#include "Simulator/DynamicBoundarySimulator.h" +#include "SPlisHSPlasH/DynamicRigidBody.h" using namespace SPH; @@ -286,31 +289,137 @@ void TimeStepDFSPH::pressureSolve() ////////////////////////////////////////////////////////////////////////// // Start solver ////////////////////////////////////////////////////////////////////////// - + Real avg_density_err = 0.0; bool chk = false; + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + Real avg_density_err_rigid = 0; + bs->computeContacts(); + bool anyRigidContact = bs->getAllContacts() != 0; + if (anyRigidContact) { + bs->computeDensityAndVolume(); + } - ////////////////////////////////////////////////////////////////////////// - // Perform solver iterations - ////////////////////////////////////////////////////////////////////////// - while ((!chk || (m_iterations < m_minIterations)) && (m_iterations < m_maxIterations)) - { - chk = true; - for (unsigned int i = 0; i < nFluids; i++) - { - FluidModel *model = sim->getFluidModel(i); - const Real density0 = model->getDensity0(); + bool rigidSolverConverged = false; + + while ((!chk || !rigidSolverConverged || (m_iterations < m_minIterations)) && (m_iterations < m_maxIterations)) { + chk = true; + rigidSolverConverged = true; + + // fluid-rigid forces + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + const Real density0 = model->getDensity0(); + const int numParticles = (int)model->numActiveParticles(); + #pragma omp for schedule(static) + for (int i = 0; i < numParticles; i++) { + computePressureAccel(fluidModelIndex, i, density0, m_simulationData.getPressureRho2Data(), true); + } + } + } + + if (sim->getDynamicBoundarySimulator() != nullptr) { + // apply force and torque for computing predicted velocity + sim->getDynamicBoundarySimulator()->updateBoundaryForces(); + + if (anyRigidContact) { + bs->computeSourceTerm(); + bs->computeDiagonalElement(); + avg_density_err_rigid = 0.0; + bs->pressureSolveIteration(avg_density_err_rigid); + if (avg_density_err_rigid > bs->getMaxDensityDeviation()) { + rigidSolverConverged = false; + } + } else { + // compute the predicted velocity and position without rigid-rigid contacts + bs->computeV_s(); + } + + // clear force and torque used to computed predicted velocity + for (int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)->getRigidBodyObject())->clearForceAndTorque(); + } + } + + + // predict density advection using the predicted velocity of fluid an rigid velocities + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int i = 0; i < (int)model->numActiveParticles(); i++) { + const Vector3r& xi = model->getPosition(i); + const Vector3r& vi = model->getVelocity(i); + const Vector3r& ai = m_simulationData.getPressureAccel(fluidModelIndex, i); + const Vector3r predictedV = vi + ai * h; + const Real& density = model->getDensity(i); + const Real density0 = model->getDensity0(); + Real densityAdv = m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i); + #ifdef USE_AVX + Scalarf8 delta_avx(0.0f); + const Vector3f8 xi_avx(xi); + Vector3f8 vi_avx(predictedV); + forall_boundary_neighbors_avx( + const Scalarf8 Vj_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVolume(0), count); + const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bs->getPredictedVelocity(pid - nFluids, 0), count); + delta_avx += Vj_avx * (vi_avx - vj_avx).dot(CubicKernel_AVX::gradW(xi_avx - xj_avx)); + ); + densityAdv += h * delta_avx.reduce(); + #else + Real delta = 0; + // density advection considering rigid bodies + forall_boundary_neighbors( + const Vector3r& vj = bm_neighbor->getVelocity(neighborIndex); + delta += bm_neighbor->getVolume(neighborIndex) * (predictedV - bs->getPredictedVelocity(pid - nFluids, neighborIndex)).dot(sim->gradW(xi - xj)); + ); + densityAdv += h * delta; + #endif + m_simulationData.getDensityAdv(fluidModelIndex, i) = densityAdv; + } + } + } + + for (unsigned int i = 0; i < nFluids; i++) { + FluidModel* model = sim->getFluidModel(i); + const Real density0 = model->getDensity0(); + + avg_density_err = 0.0; + pressureSolveIteration(i, avg_density_err); - avg_density_err = 0.0; - pressureSolveIteration(i, avg_density_err); + // Maximal allowed density fluctuation + const Real eta = m_maxError * static_cast(0.01) * density0; // maxError is given in percent - // Maximal allowed density fluctuation - const Real eta = m_maxError * static_cast(0.01) * density0; // maxError is given in percent - chk = chk && (avg_density_err <= eta); + chk = chk && (avg_density_err <= eta); + } + m_iterations++; } + } else { + ////////////////////////////////////////////////////////////////////////// + // Perform solver iterations + ////////////////////////////////////////////////////////////////////////// + while ((!chk || (m_iterations < m_minIterations)) && (m_iterations < m_maxIterations)) + { + chk = true; + for (unsigned int i = 0; i < nFluids; i++) + { + FluidModel* model = sim->getFluidModel(i); + const Real density0 = model->getDensity0(); + + avg_density_err = 0.0; + pressureSolveIteration(i, avg_density_err); + + // Maximal allowed density fluctuation + const Real eta = m_maxError * static_cast(0.01) * density0; // maxError is given in percent + chk = chk && (avg_density_err <= eta); + } - m_iterations++; + m_iterations++; + } } INCREASE_COUNTER("DFSPH - iterations", static_cast(m_iterations)); @@ -334,6 +443,13 @@ void TimeStepDFSPH::pressureSolve() } } } + + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + if (bs->getAllContacts() != 0) { + bs->applyForce(); + } + } #ifdef USE_WARMSTART for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { @@ -440,36 +556,141 @@ void TimeStepDFSPH::divergenceSolve() ////////////////////////////////////////////////////////////////////////// // Start solver ////////////////////////////////////////////////////////////////////////// - + Real avg_density_err = 0.0; bool chk = false; - ////////////////////////////////////////////////////////////////////////// - // Perform solver iterations - ////////////////////////////////////////////////////////////////////////// - while ((!chk || (m_iterationsV < 1)) && (m_iterationsV < maxIter)) - { - chk = true; - for (unsigned int i = 0; i < nFluids; i++) - { - FluidModel *model = sim->getFluidModel(i); - const Real density0 = model->getDensity0(); + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + ////////////////////////////////////////////////////////////////////////// + // Perform solver iterations (strong coupling) + ////////////////////////////////////////////////////////////////////////// + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + Real avg_density_err_rigid = 0.0; + bs->computeContacts(); + bool anyRigidContact = bs->getAllContacts() != 0; + if (anyRigidContact) { + bs->computeDensityAndVolume(); + } + bool rigidSolverConverged = false; + + while ((!chk || !rigidSolverConverged || (m_iterationsV < 1)) && (m_iterationsV < maxIter)) { + chk = true; + rigidSolverConverged = true; + // fluid-rigid forces + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + const Real density0 = model->getDensity0(); + const int numParticles = (int)model->numActiveParticles(); + #pragma omp for schedule(static) + for (int i = 0; i < numParticles; i++) { + computePressureAccel(fluidModelIndex, i, density0, m_simulationData.getPressureRho2Data(), true); + } + } + } + if (sim->getDynamicBoundarySimulator() != nullptr) { + // apply force and torque for computing predicted velocity + sim->getDynamicBoundarySimulator()->updateBoundaryForces(); + + if (anyRigidContact) { + bs->computeSourceTerm(); + bs->computeDiagonalElement(); + avg_density_err_rigid = 0.0; + bs->pressureSolveIteration(avg_density_err_rigid); + if (avg_density_err_rigid > bs->getMaxDensityDeviation()) { + rigidSolverConverged = false; + } + } else { + // compute the predicted velocity and position without rigid-rigid contacts + bs->computeV_s(); + } + // clear force and torque used to computed predicted velocity + for (int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)->getRigidBodyObject())->clearForceAndTorque(); + } + } - avg_density_err = 0.0; - divergenceSolveIteration(i, avg_density_err); + // predict density advection using the predicted velocity of fluid and rigid velocities + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int i = 0; i < (int)model->numActiveParticles(); i++) { + const Vector3r& xi = model->getPosition(i); + const Vector3r& vi = model->getVelocity(i); + const Vector3r& ai = m_simulationData.getPressureAccel(fluidModelIndex, i); + const Vector3r predictedV = vi + ai * h; + + Real densityAdv = m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i); + + #ifdef USE_AVX + Scalarf8 densityAdv_avx(0.0f); + const Vector3f8 xi_avx(xi); + Vector3f8 vi_avx(predictedV); + forall_boundary_neighbors_avx( + const Scalarf8 Vj_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVolume(0), count); + const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bs->getPredictedVelocity(pid - nFluids, 0), count); + densityAdv_avx += Vj_avx * (vi_avx - vj_avx).dot(CubicKernel_AVX::gradW(xi_avx - xj_avx)); + ); + densityAdv += densityAdv_avx.reduce(); + #else + // density advection considering rigid bodies + forall_boundary_neighbors( + densityAdv += bm_neighbor->getVolume(neighborIndex) * (predictedV - bs->getPredictedVelocity(pid - nFluids, neighborIndex)).dot(sim->gradW(xi - xj)); + ); + #endif + + m_simulationData.getDensityAdv(fluidModelIndex, i) = densityAdv; + } + } + } - // Maximal allowed density fluctuation - // use maximal density error divided by time step size - const Real eta = (static_cast(1.0) / h) * maxError * static_cast(0.01) * density0; // maxError is given in percent - chk = chk && (avg_density_err <= eta); + for (unsigned int i = 0; i < nFluids; i++) { + FluidModel* model = sim->getFluidModel(i); + const Real density0 = model->getDensity0(); + + avg_density_err = 0.0; + divergenceSolveIteration(i, avg_density_err); + + // Maximal allowed density fluctuation + // use maximal density error divided by time step size + const Real eta = (static_cast(1.0) / h) * maxError * static_cast(0.01) * density0; // maxError is given in percent + + chk = chk && (avg_density_err <= eta); + } + + m_iterationsV++; } + } else { + ////////////////////////////////////////////////////////////////////////// + // Perform solver iterations + ////////////////////////////////////////////////////////////////////////// + while ((!chk || (m_iterationsV < 1)) && (m_iterationsV < maxIter)) + { + chk = true; + for (unsigned int i = 0; i < nFluids; i++) + { + FluidModel* model = sim->getFluidModel(i); + const Real density0 = model->getDensity0(); - m_iterationsV++; - } + avg_density_err = 0.0; + divergenceSolveIteration(i, avg_density_err); - INCREASE_COUNTER("DFSPH - iterationsV", static_cast(m_iterationsV)); + // Maximal allowed density fluctuation + // use maximal density error divided by time step size + const Real eta = (static_cast(1.0) / h) * maxError * static_cast(0.01) * density0; // maxError is given in percent + chk = chk && (avg_density_err <= eta); + } + m_iterationsV++; + } + } + INCREASE_COUNTER("DFSPH - iterationsV", static_cast(m_iterationsV)); + + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { FluidModel *model = sim->getFluidModel(fluidModelIndex); @@ -491,6 +712,18 @@ void TimeStepDFSPH::divergenceSolve() } } } + + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + if (bs->getAllContacts() != 0) { + bs->applyForce(); + } + ///////////////////////////////////////////////////////////////////////// + // Update intermediate boundary velocity for strong coupling solver + ///////////////////////////////////////////////////////////////////////// + sim->getDynamicBoundarySimulator()->updateVelocities(); + } + #ifdef USE_WARMSTART_V for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { @@ -537,7 +770,7 @@ void TimeStepDFSPH::pressureSolveIteration(const unsigned int fluidModelIndex, R // (see Algorithm 3, line 7 in [BK17]) ////////////////////////////////////////////////////////////////////////// #pragma omp for schedule(static) - for (int i = 0; i < numParticles; i++) + for (int i = 0; i < numParticles; i++) { computePressureAccel(fluidModelIndex, i, density0, m_simulationData.getPressureRho2Data()); } @@ -612,10 +845,10 @@ void TimeStepDFSPH::divergenceSolveIteration(const unsigned int fluidModelIndex, { ////////////////////////////////////////////////////////////////////////// // Compute pressure accelerations using the current pressure values. - // (see Algorithm 2, line 7 in [BK17]) + // (see Algorithm 2, line 7 in [BK17]) ////////////////////////////////////////////////////////////////////////// #pragma omp for schedule(static) - for (int i = 0; i < (int)numParticles; i++) + for (int i = 0; i < (int)numParticles; i++) { computePressureAccel(fluidModelIndex, i, density0, m_simulationData.getPressureRho2VData()); } @@ -768,7 +1001,7 @@ void TimeStepDFSPH::computeDFSPHFactor(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Scalarf8 V_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVolume(0), count); @@ -870,8 +1103,12 @@ void TimeStepDFSPH::computeDensityAdv(const unsigned int fluidModelIndex, const delta += Vj * (vi - vj).dot(sim->gradW(xi - xj)); ); } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + { + m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i) = density / density0 + h * delta; + } - densityAdv = density / density0 + h*delta; + densityAdv = density / density0 + h * delta; } /** Compute rho_adv,i^(0) (see equation (9) in Section 3.2 [BK17]) @@ -933,6 +1170,10 @@ void TimeStepDFSPH::computeDensityChange(const unsigned int fluidModelIndex, con densityAdv += Vj * (vi - vj).dot(sim->gradW(xi - xj)); ); } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + { + m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i) = densityAdv; + } } /** Compute pressure accelerations using the current pressure values of the particles @@ -978,7 +1219,7 @@ void TimeStepDFSPH::computePressureAccel(const unsigned int fluidModelIndex, con ////////////////////////////////////////////////////////////////////////// if (fabs(p_rho2_i) > m_eps) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { const Scalarf8 mi_avx(model->getMass(i)); forall_boundary_neighbors_avx( @@ -1060,7 +1301,7 @@ Real TimeStepDFSPH::compute_aij_pj(const unsigned int fluidModelIndex, const uns ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Scalarf8 Vj_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVolume(0), count); @@ -1133,7 +1374,7 @@ void TimeStepDFSPH::computeDFSPHFactor(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r grad_p_j = -bm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj); @@ -1222,6 +1463,10 @@ void TimeStepDFSPH::computeDensityAdv(const unsigned int fluidModelIndex, const bm_neighbor->getPointVelocity(xj, vj); delta += Vj * (vi - vj).dot(sim->gradW(xi - xj)); ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + { + m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i) = density / density0 + h * delta; } densityAdv = density / density0 + h*delta; @@ -1278,6 +1523,10 @@ void TimeStepDFSPH::computeDensityChange(const unsigned int fluidModelIndex, con densityAdv += Vj * (vi - vj).dot(sim->gradW(xi - xj)); ); } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + { + m_simulationData.getDensityAdvNoBoundary() = densityAdv; + } } /** Compute pressure accelerations using the current pressure values of the particles @@ -1318,7 +1567,7 @@ void TimeStepDFSPH::computePressureAccel(const unsigned int fluidModelIndex, con ////////////////////////////////////////////////////////////////////////// if (fabs(p_rho2_i) > m_eps) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r grad_p_j = -bm_neighbor->getVolume(neighborIndex) * sim->gradW(xi - xj); @@ -1384,7 +1633,7 @@ Real TimeStepDFSPH::compute_aij_pj(const unsigned int fluidModelIndex, const uns ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( aij_pj += bm_neighbor->getVolume(neighborIndex) * ai.dot(sim->gradW(xi - xj)); diff --git a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.h b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.h index 07478a0e..c866806f 100644 --- a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.h +++ b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.h @@ -36,6 +36,7 @@ namespace SPH void pressureSolve(); void pressureSolveIteration(const unsigned int fluidModelIndex, Real &avg_density_err); void divergenceSolve(); + void solveRigidRigidContacts(); void divergenceSolveIteration(const unsigned int fluidModelIndex, Real &avg_density_err); void computeDensityAdv(const unsigned int fluidModelIndex, const unsigned int index, const Real h, const Real density0); void computeDensityChange(const unsigned int fluidModelIndex, const unsigned int index, const Real h); diff --git a/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp b/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp index 5d8ce94e..59d2e2a3 100644 --- a/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp +++ b/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp @@ -106,7 +106,7 @@ void DragForce_Gissler2017::step() for (unsigned int pid = 0; pid < nFluids; pid++) numNeighbors += sim->numberOfNeighbors(fluidModelIndex, pid, i); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) numNeighbors += sim->numberOfNeighbors(fluidModelIndex, pid, i); diff --git a/SPlisHSPlasH/DynamicRigidBody.h b/SPlisHSPlasH/DynamicRigidBody.h new file mode 100644 index 00000000..23531f13 --- /dev/null +++ b/SPlisHSPlasH/DynamicRigidBody.h @@ -0,0 +1,676 @@ +#ifndef __DynamicRigidBody_h__ +#define __DynamicRigidBody_h__ + +#include "Common.h" +#include "RigidBodyObject.h" +#include "TriangleMesh.h" +#include "TimeManager.h" +#include "Simulation.h" +#include "Utilities/VolumeIntegration.h" + +namespace SPH { + + class DynamicRigidBody : public RigidBodyObject { + // Some fields are from PBD::RigidBody + private: + bool m_isDynamic; + + Real m_density; + + /** mass */ + Real m_mass; + /** inverse mass */ + Real m_invMass; + /** center of mass */ + Vector3r m_x; + Vector3r m_lastX; + Vector3r m_oldX; + Vector3r m_x0; + /** center of mass velocity */ + Vector3r m_v; + Vector3r m_v0; + /** acceleration (by external forces) */ + Vector3r m_a; + /* external forces*/ + Vector3r m_force; + + /** Inertia tensor in the principal axis system: \n + * After the main axis transformation the inertia tensor is a diagonal matrix. + * So only three values are required to store the inertia tensor. These values + * are constant over time. + */ + Vector3r m_inertiaTensor; + /** 3x3 matrix, inertia tensor in world space */ + Matrix3r m_inertiaTensorW; + /** Inverse inertia tensor in body space */ + Vector3r m_inertiaTensorInverse; + /** 3x3 matrix, inverse of the inertia tensor in world space */ + Matrix3r m_inertiaTensorInverseW; + /** Quaternion that describes the rotation of the body in world space */ + Quaternionr m_q; + Quaternionr m_lastQ; + Quaternionr m_oldQ; + Quaternionr m_q0; + /** Quaternion representing the rotation of the main axis transformation + that is performed to get a diagonal inertia tensor */ + Quaternionr m_q_mat; + /** Quaternion representing the initial rotation of the geometry */ + Quaternionr m_q_initial; + /** difference of the initial translation and the translation of the main axis transformation */ + Vector3r m_x0_mat; + /** rotationMatrix = 3x3 matrix. + * Important for the transformation from world in body space and vice versa. + * When using quaternions the rotation matrix is computed out of the quaternion. + */ + Matrix3r m_rot; + /** Angular velocity, defines rotation axis and velocity (magnitude of the vector) */ + Vector3r m_omega; + Vector3r m_omega0; + /* Angular accelaration*/ + Vector3r m_a_omega; + /** external torque */ + Vector3r m_torque; + + // used to recompute mass properties, may dont need after implement the volume integrateion + Vector3r m_scale; + + // Real m_restitutionCoeff; + Real m_frictionCoeff; + + // transformation required to transform a point to local space or vice vera + Matrix3r m_transformation_R; + Vector3r m_transformation_v1; + Vector3r m_transformation_v2; + Vector3r m_transformation_R_X_v1; + protected: + + TriangleMesh m_geometry; + + public: + DynamicRigidBody() { + m_isAnimated = true; + } + + + ~DynamicRigidBody(void) + { + } + + void initBody(const Real density, const bool isDynamic, const Vector3r &position, const Quaternionr &rotation, + const Vector3r &scale, const Vector3r &velocity, const Real &friction) + { + m_isDynamic = isDynamic; + m_mass = 1.0; + m_inertiaTensor = Vector3r(1.0, 1.0, 1.0); + + m_x = position; + m_x0 = position; + m_lastX = position; + m_oldX = position; + + m_v0 = velocity; + m_v = velocity; + m_a.setZero(); + + determineMassProperties(density); + + m_q = rotation; + m_q0 = rotation; + m_lastQ = rotation; + m_oldQ = rotation; + m_rot = m_q.matrix(); + rotationUpdated(); + m_omega.setZero(); + m_omega0.setZero(); + m_torque.setZero(); + + //m_restitutionCoeff = static_cast(0.6); + m_frictionCoeff = friction; + + updateMeshTransformation(); + } + + void reset() + { + m_x = m_x0; + getOldPosition() = getPosition0(); + getLastPosition() = getPosition0(); + + m_q = m_q0; + getOldRotation() = getRotation0(); + getLastRotation() = getRotation0(); + + m_v = m_v0; + m_omega = m_omega0; + + getAcceleration().setZero(); + getForce().setZero(); + getTorque().setZero(); + + rotationUpdated(); + + updateMeshTransformation(); + + clearForceAndTorque(); + } + + void updateInverseTransformation() + { + // remove the rotation of the main axis transformation that is performed + // to get a diagonal inertia tensor since the distance function is + // evaluated in local coordinates + // + // transformation world to local: + // p_local = R_initial^T ( R_MAT R^T (p_world - x) - x_initial + x_MAT) + // + // transformation local to world: + // p_world = R R_MAT^T (R_initial p_local + x_initial - x_MAT) + x + // + m_transformation_R = (getRotationInitial().inverse() * getRotationMAT() * getRotation().inverse()).matrix(); + m_transformation_v1 = -getRotationInitial().inverse().matrix() * getPositionInitial_MAT(); + m_transformation_v2 = (getRotation()*getRotationMAT().inverse()).matrix() * getPositionInitial_MAT() + getPosition(); + m_transformation_R_X_v1 = -m_transformation_R * getPosition() + m_transformation_v1; + } + + // Determine mass and inertia tensor + void determineMassProperties(const Real density) { + m_density = density; + std::vector vd = getGeometry().getVertices(); + Utilities::VolumeIntegration vi(m_geometry.numVertices(), m_geometry.numFaces(), &m_geometry.getVertices()[0], m_geometry.getFaces().data()); + vi.compute_inertia_tensor(density); + + // Diagonalize Inertia Tensor + Eigen::SelfAdjointEigenSolver es(vi.getInertia()); + Vector3r inertiaTensor = Vector3r::Zero(); + Matrix3r R; + inertiaTensor += es.eigenvalues().x() * es.eigenvectors().col(0); + inertiaTensor += es.eigenvalues().y() * es.eigenvectors().col(1); + inertiaTensor += es.eigenvalues().z() * es.eigenvectors().col(2); + inertiaTensor.x() = std::abs(inertiaTensor.x()); + inertiaTensor.y() = std::abs(inertiaTensor.y()); + inertiaTensor.z() = std::abs(inertiaTensor.z()); + R << inertiaTensor.x(), 0, 0, + 0, inertiaTensor.y(), 0, + 0, 0, inertiaTensor.z(); + + setMass(vi.getMass()); + setInertiaTensor(inertiaTensor); + + //if (R.determinant() < 0.0) + // R = -R; + + //for (unsigned int i = 0; i < vd.size(); i++) + // vd[i] = m_rot * vd[i] + m_x0; + + //Vector3r x_MAT = vi.getCenterOfMass(); + //R = m_rot * R; + //x_MAT = m_rot * x_MAT + m_x0; + + //// rotate vertices back + //for (unsigned int i = 0; i < vd.size(); i++) + // vd[i] = R.transpose() * (vd[i] - x_MAT); + + //// set rotation + //Quaternionr qR = Quaternionr(R); + //qR.normalize(); + //m_q_mat = qR; + //m_q_initial = m_q0; + //m_x0_mat = m_x0 - x_MAT; + + //m_q0 = qR; + //m_q = m_q0; + //m_lastQ = m_q0; + //m_oldQ = m_q0; + //rotationUpdated(); + + //// set translation + //m_x0 = x_MAT; + //m_x = m_x0; + //m_lastX = m_x0; + //m_oldX = m_x0; + //updateInverseTransformation(); + } + + void rotationUpdated() + { + m_rot = m_q.matrix(); + updateInertiaW(); + updateInverseTransformation(); + } + + void updateInertiaW() + { + if (m_mass != 0.0) + { + m_inertiaTensorW = m_rot * m_inertiaTensor.asDiagonal() * m_rot.transpose(); + m_inertiaTensorInverseW = m_rot * m_inertiaTensorInverse.asDiagonal() * m_rot.transpose(); + } + } + + const Matrix3r &getTransformationR() { return m_transformation_R; } + const Vector3r &getTransformationV1() { return m_transformation_v1; } + const Vector3r &getTransformationV2() { return m_transformation_v2; } + const Vector3r &getTransformationRXV1() { return m_transformation_R_X_v1; } + + FORCE_INLINE const Vector3r& getScale() { + return m_scale; + } + + FORCE_INLINE void setMass(const Real &value) + { + m_mass = value; + if (m_mass != 0.0) + m_invMass = static_cast(1.0) / m_mass; + else + m_invMass = 0.0; + } + + FORCE_INLINE const Real &getInvMass() const + { + return m_invMass; + } + + FORCE_INLINE const Real& getDensity() const { + return m_density; + } + + FORCE_INLINE Real& getDensity() { + return m_density; + } + + FORCE_INLINE void setDensity(const Real& density) { + m_density = density; + determineMassProperties(density); + } + + FORCE_INLINE Vector3r& getPosition() { + return m_x; + } + + FORCE_INLINE Vector3r &getLastPosition() + { + return m_lastX; + } + + FORCE_INLINE const Vector3r &getLastPosition() const + { + return m_lastX; + } + + FORCE_INLINE void setLastPosition(const Vector3r &pos) + { + m_lastX = pos; + } + + FORCE_INLINE Vector3r &getOldPosition() + { + return m_oldX; + } + + FORCE_INLINE const Vector3r &getOldPosition() const + { + return m_oldX; + } + + FORCE_INLINE void setOldPosition(const Vector3r &pos) + { + m_oldX = pos; + } + + FORCE_INLINE Vector3r &getForce() + { + return m_force; + } + + FORCE_INLINE const Vector3r &getForce() const + { + return m_force; + } + + FORCE_INLINE void setForce(const Vector3r &force) + { + m_force = force; + } + + FORCE_INLINE Vector3r &getPositionInitial_MAT() + { + return m_x0_mat; + } + + FORCE_INLINE const Vector3r &getPositionInitial_MAT() const + { + return m_x0_mat; + } + + FORCE_INLINE void setPositionInitial_MAT(const Vector3r &pos) + { + m_x0_mat = pos; + } + + FORCE_INLINE Vector3r &getAcceleration() + { + return m_a; + } + + FORCE_INLINE const Vector3r &getAcceleration() const + { + return m_a; + } + + FORCE_INLINE void setAcceleration(const Vector3r &accel) + { + m_a = accel; + } + + FORCE_INLINE const Vector3r &getInertiaTensor() const + { + return m_inertiaTensor; + } + + FORCE_INLINE void setInertiaTensor(const Vector3r &value) + { + m_inertiaTensor = value; + m_inertiaTensorInverse = Vector3r(static_cast(1.0) / value[0], static_cast(1.0) / value[1], static_cast(1.0) / value[2]); + } + + FORCE_INLINE Matrix3r& getInertiaTensorW() + { + return m_inertiaTensorW; + } + + FORCE_INLINE const Matrix3r& getInertiaTensorW() const + { + return m_inertiaTensorW; + } + + FORCE_INLINE const Vector3r &getInertiaTensorInverse() const + { + return m_inertiaTensorInverse; + } + + FORCE_INLINE Matrix3r &getInertiaTensorInverseW() + { + return m_inertiaTensorInverseW; + } + + FORCE_INLINE const Matrix3r &getInertiaTensorInverseW() const + { + return m_inertiaTensorInverseW; + } + + FORCE_INLINE void setInertiaTensorInverseW(const Matrix3r &value) + { + m_inertiaTensorInverseW = value; + } + + FORCE_INLINE Quaternionr &getLastRotation() + { + return m_lastQ; + } + + FORCE_INLINE const Quaternionr &getLastRotation() const + { + return m_lastQ; + } + + FORCE_INLINE void setLastRotation(const Quaternionr &value) + { + m_lastQ = value; + } + + FORCE_INLINE Quaternionr &getOldRotation() + { + return m_oldQ; + } + + FORCE_INLINE const Quaternionr &getOldRotation() const + { + return m_oldQ; + } + + FORCE_INLINE void setOldRotation(const Quaternionr &value) + { + m_oldQ = value; + } + + FORCE_INLINE Quaternionr &getRotationMAT() + { + return m_q_mat; + } + + FORCE_INLINE const Quaternionr &getRotationMAT() const + { + return m_q_mat; + } + + FORCE_INLINE void setRotationMAT(const Quaternionr &value) + { + m_q_mat = value; + } + + FORCE_INLINE Quaternionr &getRotationInitial() + { + return m_q_initial; + } + + FORCE_INLINE const Quaternionr &getRotationInitial() const + { + return m_q_initial; + } + + FORCE_INLINE void setRotationInitial(const Quaternionr &value) + { + m_q_initial = value; + } + + FORCE_INLINE Matrix3r &getRotationMatrix() + { + return m_rot; + } + + FORCE_INLINE const Matrix3r &getRotationMatrix() const + { + return m_rot; + } + + FORCE_INLINE void setRotationMatrix(const Matrix3r &value) + { + m_rot = value; + } + + FORCE_INLINE Vector3r &getAngularAccelaration() + { + return m_a_omega; + } + + FORCE_INLINE const Vector3r &getAngularAccelaration() const + { + return m_a_omega; + } + + FORCE_INLINE void setAngularAccelaration(const Vector3r &value) + { + m_a_omega = value; + } + + FORCE_INLINE Vector3r &getAngularVelocity0() + { + return m_omega0; + } + + FORCE_INLINE const Vector3r &getAngularVelocity0() const + { + return m_omega0; + } + + FORCE_INLINE void setAngularVelocity0(const Vector3r &value) + { + m_omega0 = value; + } + + FORCE_INLINE Vector3r &getTorque() + { + return m_torque; + } + + FORCE_INLINE const Vector3r &getTorque() const + { + return m_torque; + } + + FORCE_INLINE void setTorque(const Vector3r &value) + { + m_torque = value; + } + + FORCE_INLINE Real getFrictionCoeff() const + { + return m_frictionCoeff; + } + + FORCE_INLINE void setFrictionCoeff(Real val) + { + m_frictionCoeff = val; + } + + + virtual bool isDynamic() const { + return m_isDynamic; + } + + virtual Real const getMass() const { + return m_mass; + } + virtual Vector3r const& getPosition() const { + return m_x; + } + virtual void setPosition(const Vector3r& x) { + m_x = x; + } + Vector3r const& getPosition0() const { + return m_x0; + } + void setPosition0(const Vector3r& x) { + m_x0 = x; + } + virtual Vector3r getWorldSpacePosition() const { + return m_x; + } + virtual Vector3r const& getVelocity() const { + return m_v; + } + + FORCE_INLINE Vector3r& getVelocity() { + return m_v; + } + + virtual void setVelocity(const Vector3r& v) { + if (m_isAnimated) m_v = v; + } + virtual Quaternionr const& getRotation() const { + return m_q; + } + virtual void setRotation(const Quaternionr& q) { + m_q = q; + } + Quaternionr const& getRotation0() const { + return m_q0; + } + void setRotation0(const Quaternionr& q) { + m_q0 = q; + } + virtual Matrix3r getWorldSpaceRotation() const { + return m_q.toRotationMatrix(); + } + virtual Vector3r const& getAngularVelocity() const { + return m_omega; + } + virtual void setAngularVelocity(const Vector3r& v) { + if (m_isAnimated) m_omega = v; + } + virtual void addForce(const Vector3r& f) { + m_force += f; + } + + virtual void addTorque(const Vector3r& t) { + m_torque += t; + } + + void clearForceAndTorque() { + m_force.setZero(); + m_torque.setZero(); + } + + void timeStep() { + Simulation* sim = Simulation::getCurrent(); + const Vector3r gravAccel(sim->getVecValue(Simulation::GRAVITATION)); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + // Save old values + m_lastX = m_oldX; + m_oldX = m_x; + m_lastQ = m_oldQ; + m_oldQ = m_q; + // Semi-implicit Euler + m_a = m_invMass * m_force + gravAccel; + m_v += m_a * dt; + m_x += m_v * dt; + m_a_omega = m_inertiaTensorInverseW * (m_torque - (m_omega.cross(m_inertiaTensorW * m_omega))); + m_omega += m_a_omega * dt; + Quaternionr omegaTilde(0.0, m_omega[0], m_omega[1], m_omega[2]); + m_q.coeffs() += 0.5 * (omegaTilde * m_q).coeffs() * dt; + m_q.normalize(); + rotationUpdated(); + updateMeshTransformation(); + clearForceAndTorque(); + } + + // used in DFSPH + void updateVelocity() { + Simulation* sim = Simulation::getCurrent(); + const Vector3r gravAccel(sim->getVecValue(Simulation::GRAVITATION)); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + m_a = m_invMass * m_force; + m_v += m_a * dt; + m_a_omega = m_inertiaTensorInverseW * (m_torque - (m_omega.cross(m_inertiaTensorW * m_omega))); + m_omega += m_a_omega * dt; + clearForceAndTorque(); + } + + void animate() { + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + m_x += m_v * dt; + Quaternionr angVelQ(0.0, m_omega[0], m_omega[1], m_omega[2]); + m_q.coeffs() += dt * 0.5 * (angVelQ * m_q).coeffs(); + m_q.normalize(); + updateMeshTransformation(); + } + + virtual const std::vector& getVertices() const { + return m_geometry.getVertices(); + }; + virtual const std::vector& getVertexNormals() const { + return m_geometry.getVertexNormals(); + }; + virtual const std::vector& getFaces() const { + return m_geometry.getFaces(); + }; + + void setWorldSpacePosition(const Vector3r& x) { + m_x = x; + } + void setWorldSpaceRotation(const Matrix3r& r) { + m_q = Quaternionr(r); + } + TriangleMesh& getGeometry() { + return m_geometry; + } + + virtual void updateMeshTransformation() { + m_geometry.updateMeshTransformation(m_x, m_q.toRotationMatrix()); + m_geometry.updateNormals(); + m_geometry.updateVertexNormals(); + } + }; +} + +#endif \ No newline at end of file diff --git a/SPlisHSPlasH/Emitter.cpp b/SPlisHSPlasH/Emitter.cpp index 9726e58f..131674e3 100644 --- a/SPlisHSPlasH/Emitter.cpp +++ b/SPlisHSPlasH/Emitter.cpp @@ -48,7 +48,7 @@ Vector3r Emitter::getSize(const Real width, const Real height, const int type) Vector3r size; if (type == 0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { size = { 2 * supportRadius, @@ -67,7 +67,7 @@ Vector3r Emitter::getSize(const Real width, const Real height, const int type) } else { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { // height and radius of cylinder const Real h = 2 * supportRadius; diff --git a/SPlisHSPlasH/IISPH/SimulationDataIISPH.cpp b/SPlisHSPlasH/IISPH/SimulationDataIISPH.cpp index d9c40995..e2e26ee7 100644 --- a/SPlisHSPlasH/IISPH/SimulationDataIISPH.cpp +++ b/SPlisHSPlasH/IISPH/SimulationDataIISPH.cpp @@ -25,6 +25,7 @@ void SimulationDataIISPH::init() m_pressure.resize(nModels); m_lastPressure.resize(nModels); m_pressureAccel.resize(nModels); + m_density_adv_no_boundary.resize(nModels); for (unsigned int i = 0; i < nModels; i++) { FluidModel *fm = sim->getFluidModel(i); @@ -35,6 +36,7 @@ void SimulationDataIISPH::init() m_pressure[i].resize(fm->numParticles(), 0.0); m_lastPressure[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); } } @@ -52,6 +54,7 @@ void SimulationDataIISPH::cleanup() m_pressure[i].clear(); m_lastPressure[i].clear(); m_pressureAccel[i].clear(); + m_density_adv_no_boundary[i].clear(); } m_aii.clear(); m_dii.clear(); @@ -60,6 +63,7 @@ void SimulationDataIISPH::cleanup() m_pressure.clear(); m_lastPressure.clear(); m_pressureAccel.clear(); + m_density_adv_no_boundary.clear(); } void SimulationDataIISPH::reset() @@ -96,6 +100,7 @@ void SimulationDataIISPH::performNeighborhoodSearchSort() d.sort_field(&m_pressure[i][0]); d.sort_field(&m_lastPressure[i][0]); d.sort_field(&m_pressureAccel[i][0]); + d.sort_field(&m_density_adv_no_boundary[i][0]); } } } diff --git a/SPlisHSPlasH/IISPH/SimulationDataIISPH.h b/SPlisHSPlasH/IISPH/SimulationDataIISPH.h index f3aaef9c..636b1256 100644 --- a/SPlisHSPlasH/IISPH/SimulationDataIISPH.h +++ b/SPlisHSPlasH/IISPH/SimulationDataIISPH.h @@ -24,6 +24,8 @@ namespace SPH std::vector> m_dii; std::vector> m_dij_pj; std::vector> m_density_adv; + // for Gissler2019 strong coupling + std::vector> m_density_adv_no_boundary; std::vector> m_pressure; std::vector> m_lastPressure; std::vector> m_pressureAccel; @@ -108,6 +110,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 getPressure(const unsigned int fluidIndex, const unsigned int i) const { return m_pressure[fluidIndex][i]; diff --git a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp index 55bf8d2c..124807f1 100644 --- a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp +++ b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp @@ -9,7 +9,9 @@ #include "SPlisHSPlasH/BoundaryModel_Akinci2012.h" #include "SPlisHSPlasH/BoundaryModel_Koschier2017.h" #include "SPlisHSPlasH/BoundaryModel_Bender2019.h" - +#include "SPlisHSPlasH/StrongCouplingBoundarySolver.h" +#include "Simulator/DynamicBoundarySimulator.h" +#include "SPlisHSPlasH/DynamicRigidBody.h" using namespace SPH; using namespace std; @@ -87,12 +89,23 @@ void TimeStepIISPH::step() STOP_TIMING_AVG; START_TIMING("pressureSolve"); - pressureSolve(); + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + pressureSolveStrongCoupling(); + } else { + pressureSolve(); + } STOP_TIMING_AVG; for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++) integration(fluidModelIndex); + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + if (bs->getAllContacts() != 0) { + bs->applyForce(); + } + } + sim->emitParticles(); sim->animateParticles(); @@ -151,7 +164,7 @@ void TimeStepIISPH::predictAdvection(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( dii -= bm_neighbor->getVolume(neighborIndex) / density2 * sim->gradW(xi - xj); @@ -202,7 +215,7 @@ void TimeStepIISPH::predictAdvection(const unsigned int fluidModelIndex) densityAdv += h*bm_neighbor->getVolume(neighborIndex) * (vi - vj).dot(sim->gradW(xi - xj)); ); } - else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) { forall_density_maps( Vector3r vj; @@ -217,6 +230,11 @@ void TimeStepIISPH::predictAdvection(const unsigned int fluidModelIndex) bm_neighbor->getPointVelocity(xj, vj); densityAdv += h*Vj * (vi - vj).dot(sim->gradW(xi - xj)); ); + } + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + { + // for strong coupling, the density advection due to boundary particles are computed in pressure solve iterations. + m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i) = densityAdv; } const Real &pressure = m_simulationData.getPressure(fluidModelIndex, i); @@ -245,7 +263,7 @@ void TimeStepIISPH::predictAdvection(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r kernel = sim->gradW(xi - xj); @@ -302,6 +320,126 @@ void TimeStepIISPH::pressureSolve() INCREASE_COUNTER("IISPH - iterations", static_cast(m_iterations)); } +void TimeStepIISPH::pressureSolveStrongCoupling() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + Real avg_density_err = 0; + Real avg_density_err_rigid = 0; + m_iterations = 0; + bool chk = false; + + bs->computeContacts(); + bool anyRigidContact = bs->getAllContacts() != 0; + if (anyRigidContact) { + bs->computeDensityAndVolume(); + } + bool rigidSolverConverged = false; + + while ((!chk || !rigidSolverConverged || (m_iterations < m_minIterations)) && (m_iterations < m_maxIterations)) { + chk = true; + rigidSolverConverged = true; + + // fluid-rigid forces + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int i = 0; i < (int)model->numActiveParticles(); i++) { + const Real density0 = model->getDensity0(); + const Real density_i = model->getDensity(i) / density0; + const Real density2 = density_i * density_i; + const Real dpi = m_simulationData.getPressure(fluidModelIndex, i) / density2; + const Vector3r& xi = model->getPosition(i); + const Vector3r& vi = model->getVelocity(i); + + Vector3r& ai = m_simulationData.getPressureAccel(fluidModelIndex, i); + ai.setZero(); + + // fluid pressure + forall_fluid_neighbors( + const Real density_j = fm_neighbor->getDensity(neighborIndex) / fm_neighbor->getDensity0(); + const Real densityj2 = density_j * density_j; + const Real dpj = m_simulationData.getPressure(pid, neighborIndex) / densityj2; + ai -= fm_neighbor->getVolume(neighborIndex) * (dpi + fm_neighbor->getDensity0() / density0 * dpj) * sim->gradW(xi - xj); + ); + + // fluid-rigid + forall_boundary_neighbors( + const Vector3r a = bm_neighbor->getVolume(neighborIndex) * dpi * sim->gradW(xi - xj); + ai -= a; + bm_neighbor->addForce(xj, model->getMass(i) * a); + ); + } + } + } + + // apply force and torque for computing predicted velocity + if (sim->getDynamicBoundarySimulator() != nullptr) { + sim->getDynamicBoundarySimulator()->updateBoundaryForces(); + + if (anyRigidContact) { + bs->computeSourceTerm(); + bs->computeDiagonalElement(); + avg_density_err_rigid = 0.0; + bs->pressureSolveIteration(avg_density_err_rigid); + if (avg_density_err_rigid > bs->getMaxDensityDeviation()) { + rigidSolverConverged = false; + } + } else { + // compute the predicted velocity and position without rigid-rigid contacts + bs->computeV_s(); + } + + // clear force and torque used to computed predicted velocity + for (int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)->getRigidBodyObject())->clearForceAndTorque(); + } + } + + // predict density advection using the predicted velocity of fluid and rigid velocities + for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++) { + FluidModel* model = sim->getFluidModel(fluidModelIndex); + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int i = 0; i < (int)model->numActiveParticles(); i++) { + const Vector3r& xi = model->getPosition(i); + const Vector3r& vi = model->getVelocity(i); + const Vector3r& ai = m_simulationData.getPressureAccel(fluidModelIndex, i); + const Vector3r predictedV = vi + ai * dt; + + Real densityAdv = m_simulationData.getDensityAdvNoBoundary(fluidModelIndex, i); + + // density advection considering rigid bodies + forall_boundary_neighbors( + densityAdv += dt * bm_neighbor->getVolume(neighborIndex) * (predictedV - bs->getPredictedVelocity(pid - nFluids, neighborIndex)).dot(sim->gradW(xi - xj)); + ); + m_simulationData.getDensityAdv(fluidModelIndex, i) = densityAdv; + } + } + } + + for (unsigned int i = 0; i < nFluids; i++) { + FluidModel* model = sim->getFluidModel(i); + const Real density0 = model->getDensity0(); + + avg_density_err = 0.0; + pressureSolveIteration(i, avg_density_err); + + // Maximal allowed density fluctuation + const Real eta = m_maxError * static_cast(0.01) * density0; // maxError is given in percent + chk = chk && (avg_density_err <= eta); + } + + m_iterations++; + } + + INCREASE_COUNTER("IISPH - iterations", static_cast(m_iterations)); +} + void TimeStepIISPH::pressureSolveIteration(const unsigned int fluidModelIndex, Real &avg_density_err) { Simulation *sim = Simulation::getCurrent(); @@ -377,7 +515,7 @@ void TimeStepIISPH::pressureSolveIteration(const unsigned int fluidModelIndex, R ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( sum += bm_neighbor->getVolume(neighborIndex) * m_simulationData.getDij_pj(fluidModelIndex, i).dot(sim->gradW(xi - xj)); @@ -492,7 +630,7 @@ void TimeStepIISPH::computePressureAccels(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r a = bm_neighbor->getVolume(neighborIndex) * dpi * sim->gradW(xi - xj); diff --git a/SPlisHSPlasH/IISPH/TimeStepIISPH.h b/SPlisHSPlasH/IISPH/TimeStepIISPH.h index 77b40292..43779e86 100644 --- a/SPlisHSPlasH/IISPH/TimeStepIISPH.h +++ b/SPlisHSPlasH/IISPH/TimeStepIISPH.h @@ -24,6 +24,7 @@ namespace SPH void predictAdvection(const unsigned int fluidModelIndex); void pressureSolve(); + void pressureSolveStrongCoupling(); void pressureSolveIteration(const unsigned int fluidModelIndex, Real &avg_density_err); void integration(const unsigned int fluidModelIndex); diff --git a/SPlisHSPlasH/PBF/TimeStepPBF.cpp b/SPlisHSPlasH/PBF/TimeStepPBF.cpp index eafb0553..860f5a9e 100644 --- a/SPlisHSPlasH/PBF/TimeStepPBF.cpp +++ b/SPlisHSPlasH/PBF/TimeStepPBF.cpp @@ -251,7 +251,7 @@ void TimeStepPBF::pressureSolveIteration(const unsigned int fluidModelIndex, Rea density += rho; ); } - else // Bender2019 + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) { forall_volume_maps( density += Vj * sim->W(xi - xj); @@ -297,7 +297,7 @@ void TimeStepPBF::pressureSolveIteration(const unsigned int fluidModelIndex, Rea gradC_i -= gradRho; ); } - else // Bender2019 + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) { forall_volume_maps( const Vector3r gradC_j = -Vj * sim->gradW(xi - xj); @@ -360,7 +360,7 @@ void TimeStepPBF::pressureSolveIteration(const unsigned int fluidModelIndex, Rea bm_neighbor->addForce(xj, model->getMass(i) * dx * invH2); ); } - else // Bender2019 + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) { forall_volume_maps( const Vector3r gradC_j = -Vj * sim->gradW(xi - xj); diff --git a/SPlisHSPlasH/Simulation.cpp b/SPlisHSPlasH/Simulation.cpp index 79797b87..5994eca1 100644 --- a/SPlisHSPlasH/Simulation.cpp +++ b/SPlisHSPlasH/Simulation.cpp @@ -13,6 +13,7 @@ #include "BoundaryModel_Akinci2012.h" #include "BoundaryModel_Bender2019.h" #include "BoundaryModel_Koschier2017.h" +#include "StrongCouplingBoundarySolver.h" using namespace SPH; @@ -59,6 +60,7 @@ int Simulation::BOUNDARY_HANDLING_METHOD = -1; int Simulation::ENUM_AKINCI2012 = -1; int Simulation::ENUM_KOSCHIER2017 = -1; int Simulation::ENUM_BENDER2019 = -1; +int Simulation::ENUM_GISSLER2019 = -1; Simulation::Simulation () @@ -83,6 +85,7 @@ Simulation::Simulation () m_animationFieldSystem = new AnimationFieldSystem(); m_boundaryHandlingMethod = static_cast(BoundaryHandlingMethods::Bender2019); + m_dynamicBoundarySimulator = nullptr; } Simulation::~Simulation () @@ -265,6 +268,7 @@ void Simulation::initParameters() enumParam->addEnumValue("Akinci et al. 2012", ENUM_AKINCI2012); enumParam->addEnumValue("Koschier and Bender 2017", ENUM_KOSCHIER2017); enumParam->addEnumValue("Bender et al. 2019", ENUM_BENDER2019); + enumParam->addEnumValue("Gissler et al. 2019", ENUM_GISSLER2019); enumParam->setReadOnly(true); } @@ -380,7 +384,7 @@ void Simulation::setKernel(int val) m_kernelFct = WendlandQuinticC2Kernel2D::W; } } - if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) updateBoundaryVolume(); } @@ -429,7 +433,7 @@ void Simulation::updateTimeStepSizeCFL() } // boundary particles - if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { for (unsigned int i = 0; i < numberOfBoundaryModels(); i++) { @@ -508,8 +512,14 @@ void Simulation::reset() for (unsigned int i = 0; i < numberOfBoundaryModels(); i++) getBoundaryModel(i)->reset(); - if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { updateBoundaryVolume(); + } + + if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver::getCurrent()->reset(); + } + if (m_timeStep) m_timeStep->reset(); @@ -616,6 +626,9 @@ void Simulation::performNeighborhoodSearchSort() BoundaryModel *bm = getBoundaryModel(i); bm->performNeighborhoodSearchSort(); } + if (getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + StrongCouplingBoundarySolver::getCurrent()->performNeighborhoodSearchSort(); + } #ifdef USE_DEBUG_TOOLS m_debugTools->performNeighborhoodSearchSort(); #endif @@ -733,6 +746,9 @@ void Simulation::updateBoundaryVolume() for (unsigned int j = numberOfFluidModels(); j < m_neighborhoodSearch->point_sets().size(); j++) m_neighborhoodSearch->set_active(i, j, true); } + + // Activate all point sets for strong fluid-rigid coupling + m_neighborhoodSearch->set_active(true); } void SPH::Simulation::saveState(BinaryFileWriter &binWriter) diff --git a/SPlisHSPlasH/Simulation.h b/SPlisHSPlasH/Simulation.h index d0f0260c..67ddcff4 100644 --- a/SPlisHSPlasH/Simulation.h +++ b/SPlisHSPlasH/Simulation.h @@ -9,6 +9,7 @@ #include "BoundaryModel.h" #include "AnimationFieldSystem.h" #include "Utilities/FileSystem.h" + #ifdef USE_DEBUG_TOOLS #include "SPlisHSPlasH/Utilities/DebugTools.h" #endif @@ -167,8 +168,9 @@ for (unsigned int pid = 0; pid < nBoundaries; pid++) \ namespace SPH { + class DynamicBoundarySimulator; enum class SimulationMethods { WCSPH = 0, PCISPH, PBF, IISPH, DFSPH, PF, ICSPH, NumSimulationMethods }; - enum class BoundaryHandlingMethods { Akinci2012 = 0, Koschier2017, Bender2019, NumSimulationMethods }; + enum class BoundaryHandlingMethods { Akinci2012 = 0, Koschier2017, Bender2019, Gissler2019, NumSimulationMethods }; /** \brief Class to manage the current simulation time and the time step size. * This class is a singleton. @@ -220,6 +222,7 @@ namespace SPH static int ENUM_AKINCI2012; static int ENUM_KOSCHIER2017; static int ENUM_BENDER2019; + static int ENUM_GISSLER2019; typedef PrecomputedKernel PrecomputedCubicKernel; @@ -306,6 +309,7 @@ namespace SPH bool m_enableZSort; std::function m_simulationMethodChanged; int m_boundaryHandlingMethod; + DynamicBoundarySimulator* m_dynamicBoundarySimulator; std::string m_cachePath; bool m_useCache; std::vector m_dragMethods; @@ -366,6 +370,9 @@ namespace SPH BoundaryHandlingMethods getBoundaryHandlingMethod() const { return (BoundaryHandlingMethods) m_boundaryHandlingMethod; } void setBoundaryHandlingMethod(BoundaryHandlingMethods val) { m_boundaryHandlingMethod = (int) val; } + DynamicBoundarySimulator* getDynamicBoundarySimulator() const { return m_dynamicBoundarySimulator; } + void setDynamicBoundarySimulator(DynamicBoundarySimulator* simulator) { m_dynamicBoundarySimulator = simulator; } + int getKernel() const { return m_kernelMethod; } void setKernel(int val); int getGradKernel() const { return m_gradKernelMethod; } diff --git a/SPlisHSPlasH/StrongCouplingBoundarySolver.cpp b/SPlisHSPlasH/StrongCouplingBoundarySolver.cpp new file mode 100644 index 00000000..44ea7b06 --- /dev/null +++ b/SPlisHSPlasH/StrongCouplingBoundarySolver.cpp @@ -0,0 +1,753 @@ +#include "SPlisHSPlasH/BoundaryModel_Akinci2012.h" +#include "StrongCouplingBoundarySolver.h" +#include "TimeManager.h" +#include "DynamicRigidBody.h" +#include "Simulator/DynamicBoundarySimulator.h" +#include "SPlisHSPlasH/SPHKernels.h" + +using namespace SPH; +using namespace GenParam; + +int StrongCouplingBoundarySolver::KERNEL_METHOD = -1; +int StrongCouplingBoundarySolver::GRAD_KERNEL_METHOD = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_CUBIC = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_WENDLANDQUINTICC2 = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_POLY6 = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_SPIKY = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_CUBIC_2D = -1; +int StrongCouplingBoundarySolver::ENUM_KERNEL_WENDLANDQUINTICC2_2D = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_CUBIC = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_WENDLANDQUINTICC2 = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_POLY6 = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_SPIKY = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_CUBIC_2D = -1; +int StrongCouplingBoundarySolver::ENUM_GRADKERNEL_WENDLANDQUINTICC2_2D = -1; +StrongCouplingBoundarySolver* StrongCouplingBoundarySolver::current = nullptr; + +StrongCouplingBoundarySolver::StrongCouplingBoundarySolver() : + m_density(), + m_pressureGrad(), + m_v_s(), + m_s(), + m_pressure(), + m_v_rr(), + m_minus_rho_div_v_s(), + m_minus_rho_div_v_rr(), + m_diagonalElement(), + m_artificialVolume(), + m_predictedVelocity(), + m_predictedPosition(), + m_v_rr_body(), + m_omega_rr_body(), + m_bodyContacts() +{ + m_maxIterations = 100; + m_minIterations = 2; + m_maxDensityDeviation = 0.001; + + m_kernelMethod = 0; + m_gradKernelMethod = 0; + m_kernelFct = CubicKernel::W; + m_gradKernelFct = CubicKernel::gradW; + m_W_zero = CubicKernel::W_zero(); + + initParameters(); + + Simulation* sim = Simulation::getCurrent(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + + // export fields should be added in resize() + resize(nBoundaries); +} + + +void SPH::StrongCouplingBoundarySolver::resize(unsigned int size) { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + + m_restDensity = 1000; + m_density.resize(nBoundaries); + m_v_s.resize(nBoundaries); + m_s.resize(nBoundaries); + m_pressure.resize(nBoundaries); + m_v_rr.resize(nBoundaries); + m_minus_rho_div_v_s.resize(nBoundaries); + m_minus_rho_div_v_rr.resize(nBoundaries); + m_diagonalElement.resize(nBoundaries); + m_pressureGrad.resize(nBoundaries); + m_artificialVolume.resize(nBoundaries); + m_predictedVelocity.resize(nBoundaries); + m_predictedPosition.resize(nBoundaries); + m_v_rr_body.resize(nBoundaries, Vector3r::Zero()); + m_omega_rr_body.resize(nBoundaries, Vector3r::Zero()); + m_bodyContacts.resize(nBoundaries, 0); + + for (unsigned int i = 0; i < nBoundaries; i++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModel(i)); + m_density[i].resize(bm->numberOfParticles(), m_restDensity); + m_v_s[i].resize(bm->numberOfParticles(), Vector3r::Zero()); + m_s[i].resize(bm->numberOfParticles(), 0.0); + m_pressure[i].resize(bm->numberOfParticles(), 0.0); + m_v_rr[i].resize(bm->numberOfParticles(), Vector3r::Zero()); + m_minus_rho_div_v_s[i].resize(bm->numberOfParticles(), 0.0); + m_minus_rho_div_v_rr[i].resize(bm->numberOfParticles(), 0.0); + m_diagonalElement[i].resize(bm->numberOfParticles(), 0.0); + m_pressureGrad[i].resize(bm->numberOfParticles(), Vector3r::Zero()); + m_artificialVolume[i].resize(bm->numberOfParticles(), 0.0); + m_predictedVelocity[i].resize(bm->numberOfParticles(), Vector3r::Zero()); + m_predictedPosition[i].resize(bm->numberOfParticles(), Vector3r::Zero()); + + //bm->addField({ "particle position", FieldType::Vector3, [sim, bm, i](const unsigned int j)->Vector3r* {return &bm->getPosition(j); } }); + bm->addField({ "body position", FieldType::Vector3, [bm](const unsigned int j)->Vector3r* {return &static_cast(bm->getRigidBodyObject())->getPosition(); } }); + } +} + + +void SPH::StrongCouplingBoundarySolver::reset() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nBoundaries = sim->numberOfBoundaryModels(); + + for (unsigned int i = 0; i < nBoundaries; i++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModel(i)); + for (int j = 0; j < bm->numberOfParticles(); j++) { + m_density[i][j] = m_restDensity; + m_v_s[i][j] = Vector3r::Zero(); + m_s[i][j] = 0.0; + m_pressure[i][j] = 0.0; + m_v_rr[i][j] = Vector3r::Zero(); + m_minus_rho_div_v_s[i][j] = 0.0; + m_minus_rho_div_v_rr[i][j] = 0.0; + m_diagonalElement[i][j] = 0.0; + m_pressureGrad[i][j] = Vector3r::Zero(); + m_artificialVolume[i][j] = 0.0; + m_predictedVelocity[i][j] = Vector3r::Zero(); + m_predictedPosition[i][j] = Vector3r::Zero(); + } + m_v_rr_body[i] = Vector3r::Zero(); + m_omega_rr_body[i] = Vector3r::Zero(); + m_bodyContacts[i] = 0; + } +} + +StrongCouplingBoundarySolver::~StrongCouplingBoundarySolver() { + current = nullptr; + for (unsigned int i = 0; i < Simulation::getCurrent()->numberOfBoundaryModels(); i++) { + m_density[i].clear(); + m_v_s[i].clear(); + m_s[i].clear(); + m_pressure[i].clear(); + m_v_rr[i].clear(); + m_minus_rho_div_v_s[i].clear(); + m_minus_rho_div_v_rr[i].clear(); + m_diagonalElement[i].clear(); + m_pressureGrad[i].clear(); + m_artificialVolume[i].clear(); + m_predictedVelocity[i].clear(); + m_predictedPosition[i].clear(); + } + m_density.clear(); + m_v_s.clear(); + m_s.clear(); + m_pressure.clear(); + m_v_rr.clear(); + m_minus_rho_div_v_s.clear(); + m_minus_rho_div_v_rr.clear(); + m_diagonalElement.clear(); + m_pressureGrad.clear(); + m_artificialVolume.clear(); + m_predictedVelocity.clear(); + m_predictedPosition.clear(); + m_v_rr_body.clear(); + m_omega_rr_body.clear(); + m_bodyContacts.clear(); +} + +StrongCouplingBoundarySolver* StrongCouplingBoundarySolver::getCurrent() { + if (current == nullptr) { + current = new StrongCouplingBoundarySolver(); + } else { + if (Simulation::getCurrent()->numberOfBoundaryModels() != current->m_density.size()) { + current->resize(Simulation::getCurrent()->numberOfBoundaryModels()); + } + } + return current; +} + +void StrongCouplingBoundarySolver::initParameters() { + ParameterObject::initParameters(); + + Simulation* sim = Simulation::getCurrent(); + ParameterBase::GetFunc getKernelFct = std::bind(&StrongCouplingBoundarySolver::getKernel, this); + ParameterBase::SetFunc setKernelFct = std::bind(&StrongCouplingBoundarySolver::setKernel, this, std::placeholders::_1); + KERNEL_METHOD = createEnumParameter("kernel", "Kernel", getKernelFct, setKernelFct); + setGroup(KERNEL_METHOD, "Boundary Model|Kernel"); + setDescription(KERNEL_METHOD, "Kernel function used in the boundary solver."); + EnumParameter* enumParam = static_cast(getParameter(KERNEL_METHOD)); + if (!sim->is2DSimulation()) { + enumParam->addEnumValue("Cubic spline", ENUM_KERNEL_CUBIC); + enumParam->addEnumValue("Wendland quintic C2", ENUM_KERNEL_WENDLANDQUINTICC2); + enumParam->addEnumValue("Poly6", ENUM_KERNEL_POLY6); + enumParam->addEnumValue("Spiky", ENUM_KERNEL_SPIKY); + } else { + enumParam->addEnumValue("Cubic spline 2D", ENUM_KERNEL_CUBIC_2D); + enumParam->addEnumValue("Wendland quintic C2 2D", ENUM_KERNEL_WENDLANDQUINTICC2_2D); + } + + ParameterBase::GetFunc getGradKernelFct = std::bind(&StrongCouplingBoundarySolver::getGradKernel, this); + ParameterBase::SetFunc setGradKernelFct = std::bind(&StrongCouplingBoundarySolver::setGradKernel, this, std::placeholders::_1); + GRAD_KERNEL_METHOD = createEnumParameter("gradKernel", "Gradient of kernel", getGradKernelFct, setGradKernelFct); + setGroup(GRAD_KERNEL_METHOD, "Boundary Model|Kernel"); + setDescription(GRAD_KERNEL_METHOD, "Gradient of the kernel function used in the boundary solver."); + enumParam = static_cast(getParameter(GRAD_KERNEL_METHOD)); + if (!sim->is2DSimulation()) { + enumParam->addEnumValue("Cubic spline", ENUM_GRADKERNEL_CUBIC); + enumParam->addEnumValue("Wendland quintic C2", ENUM_GRADKERNEL_WENDLANDQUINTICC2); + enumParam->addEnumValue("Poly6", ENUM_GRADKERNEL_POLY6); + enumParam->addEnumValue("Spiky", ENUM_GRADKERNEL_SPIKY); + } else { + enumParam->addEnumValue("Cubic spline 2D", ENUM_GRADKERNEL_CUBIC_2D); + enumParam->addEnumValue("Wendland quintic C2 2D", ENUM_GRADKERNEL_WENDLANDQUINTICC2_2D); + } +} + +void StrongCouplingBoundarySolver::performNeighborhoodSearchSort() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nModels = sim->numberOfBoundaryModels(); + + for (unsigned int i = 0; i < nModels; i++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModel(i)); + + const unsigned int numPart = bm->numberOfParticles(); + if (numPart != 0) { + auto const& d = sim->getNeighborhoodSearch()->point_set(bm->getPointSetIndex()); + d.sort_field(&m_density[i][0]); + d.sort_field(&m_v_s[i][0]); + d.sort_field(&m_s[i][0]); + d.sort_field(&m_pressure[i][0]); + d.sort_field(&m_v_rr[i][0]); + d.sort_field(&m_minus_rho_div_v_rr[i][0]); + d.sort_field(&m_diagonalElement[i][0]); + d.sort_field(&m_pressureGrad[i][0]); + d.sort_field(&m_artificialVolume[i][0]); + d.sort_field(&m_predictedVelocity[i][0]); + d.sort_field(&m_predictedPosition[i][0]); + } + } +} + +void StrongCouplingBoundarySolver::computeContacts() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + m_contactsAllBodies = 0; + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + unsigned int contacts = 0; + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + int bmIndex = boundaryPointSetIndex - nFluids; + if (bm->getRigidBodyObject()->isDynamic()) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) reduction(+:contacts) + for (int r = 0; r < bm->numberOfParticles(); r++) { + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + if (pid != boundaryPointSetIndex) { + if (sim->numberOfNeighbors(boundaryPointSetIndex, pid, r) > 0) { + contacts++; + break; + } + } + } + } + } + } + m_bodyContacts[bmIndex] = contacts; + m_contactsAllBodies += contacts; + } +} + + +void StrongCouplingBoundarySolver::computeDensityAndVolume() { + Simulation* sim = Simulation::getCurrent(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + + // Compute density and artificial volume for all particles + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + int bmIndex = boundaryPointSetIndex - nFluids; + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm->numberOfParticles(); r++) { + const Real gamma = static_cast(1); + if (bm->getRigidBodyObject()->isDynamic() && m_bodyContacts[bmIndex] > 0) { + // compute density for particle r + Real particleDensity = getRestDensity() * bm->getVolume(r) * W_zero(); + // iterate over all rigid bodies + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + BoundaryModel_Akinci2012* bm_neighbor = static_cast(sim->getBoundaryModelFromPointSet(pid)); + for (unsigned int j = 0; j < sim->numberOfNeighbors(boundaryPointSetIndex, pid, r); j++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, pid, r, j); + particleDensity += getRestDensity() * bm->getVolume(r) * W(bm->getPosition(r) - bm_neighbor->getPosition(k)); + } + } + setDensity(bmIndex, r, std::max(particleDensity, getRestDensity())); + if (getDensity(bmIndex, r) > getRestDensity()) { + setArtificialVolume(bmIndex, r, gamma * getRestDensity() * bm->getVolume(r) / getDensity(bmIndex, r)); + } else { + setArtificialVolume(bmIndex, r, gamma * bm->getVolume(r)); + } + } else { + setDensity(bmIndex, r, getRestDensity()); + setArtificialVolume(bmIndex, r, gamma * bm->getVolume(r)); + } + } + } + } +} + +void StrongCouplingBoundarySolver::computeV_s() { + Simulation* sim = Simulation::getCurrent(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + + // Compute v_s for all particles + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + DynamicRigidBody* rb = dynamic_cast(bm->getRigidBodyObject()); + if (rb != nullptr && rb->isDynamic() && m_bodyContacts[boundaryPointSetIndex - nFluids] > 0) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm->numberOfParticles(); r++) { + Vector3r gravForce = rb->getMass() * Vector3r(sim->getVecValue(Simulation::GRAVITATION)); + // the rb->getForce() contains only fluid-rigid forces + Vector3r F_R = rb->getForce() + gravForce; + Vector3r pos = bm->getPosition(r) - rb->getPosition(); + // compute v_s + setV_s(boundaryPointSetIndex - nFluids, r, rb->getVelocity() + dt * rb->getInvMass() * F_R + + (rb->getAngularVelocity() + rb->getInertiaTensorInverseW() * (dt * rb->getTorque() + dt * (rb->getInertiaTensorW() * rb->getAngularVelocity()).cross(rb->getAngularVelocity()))).cross(pos)); + // compute the predicted velocity first here, since for cases without rigid-rigid contacts, the computeSourceTermRHS will not be called and V_rr is 0 + // but we still need the predicted velocity for fluid-rigid strong coupling + setPredictedVelocity(boundaryPointSetIndex - nFluids, r, getV_s(boundaryPointSetIndex - nFluids, r)); + setPredictedPosition(boundaryPointSetIndex - nFluids, r, bm->getPosition(r) + (bm->getVelocity(r) + getV_s(boundaryPointSetIndex - nFluids, r)) / 2 * dt); + } + } + } + } +} + +void StrongCouplingBoundarySolver::computeSourceTerm() { + Simulation* sim = Simulation::getCurrent(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + + computeV_s(); + + // compute source term s for all particles + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + int bmIndex = boundaryPointSetIndex - nFluids; + if (bm->getRigidBodyObject()->isDynamic() && m_bodyContacts[bmIndex] > 0) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm->numberOfParticles(); r++) { + Real rho_div_v_s = 0; + // iterate over all rigid bodies except bm, since the divergence of particles in the same rigid body should be 0. + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + if (boundaryPointSetIndex != pid) { + BoundaryModel_Akinci2012* bm_neighbor = static_cast(sim->getBoundaryModelFromPointSet(pid)); + int neighborIndex = pid - nFluids; + for (unsigned int j = 0; j < sim->numberOfNeighbors(boundaryPointSetIndex, pid, r); j++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, pid, r, j); + // sum up divergence of v_s + rho_div_v_s += getArtificialVolume(neighborIndex, k) * getDensity(neighborIndex, k) * (getV_s(neighborIndex, k) - getV_s(bmIndex, r)).dot(gradW(bm->getPosition(r) - bm_neighbor->getPosition(k))); + } + } + } + Real s = (getRestDensity() - getDensity(bmIndex, r)) / dt + rho_div_v_s; + setSourceTerm(bmIndex, r, s); + // used to computed the density deviation + setMinus_rho_div_v_s(bmIndex, r, -rho_div_v_s); + } + } + } + } +} + +void StrongCouplingBoundarySolver::computeDiagonalElement() { + Simulation* sim = Simulation::getCurrent(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + for (unsigned int pointSetIndex_r = nFluids; pointSetIndex_r < sim->numberOfPointSets(); pointSetIndex_r++) { + BoundaryModel_Akinci2012* bm_r = static_cast(sim->getBoundaryModelFromPointSet(pointSetIndex_r)); + int index_r = pointSetIndex_r - nFluids; + DynamicRigidBody* rb = dynamic_cast(bm_r->getRigidBodyObject()); + if (rb != nullptr && rb->isDynamic()) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm_r->numberOfParticles(); r++) { + if (m_bodyContacts[index_r] == 0) { + setDiagonalElement(index_r, r, 0); + } else { + Vector3r pos_r = bm_r->getPosition(r) - rb->getPosition(); + Vector3r grad_p_b = Vector3r::Zero(); + const Real density_r2 = getDensity(index_r, r) * getDensity(index_r, r); + // rk are neighboring rigid particles of r of other rigid bodies k + for (unsigned int pointSetIndex_rk = nFluids; pointSetIndex_rk < sim->numberOfPointSets(); pointSetIndex_rk++) { + if (pointSetIndex_rk != pointSetIndex_r) { + BoundaryModel_Akinci2012* bm_rk = static_cast(sim->getBoundaryModelFromPointSet(pointSetIndex_rk)); + int index_rk = pointSetIndex_rk - nFluids; + DynamicRigidBody* rb_rk = static_cast(bm_rk->getRigidBodyObject()); + for (unsigned int n = 0; n < sim->numberOfNeighbors(pointSetIndex_r, pointSetIndex_rk, r); n++) { + const unsigned int rk = sim->getNeighbor(pointSetIndex_r, pointSetIndex_rk, r, n); + grad_p_b += getArtificialVolume(index_rk, rk) * getDensity(index_rk, rk) / density_r2 * gradW(bm_r->getPosition(r) - bm_rk->getPosition(rk)); + } + } + } + grad_p_b *= getDensity(index_r, r); + Vector3r v_b = -dt * rb->getInvMass() * getArtificialVolume(index_r, r) * grad_p_b; + Vector3r omega_b = -dt * rb->getInertiaTensorInverseW() * getArtificialVolume(index_r, r) * pos_r.cross(grad_p_b); + Vector3r v_b_r = v_b + omega_b.cross(pos_r); + Real b_r = 0; + // for all neighboring rigid bodies + for (unsigned int pointSetIndex_rk = nFluids; pointSetIndex_rk < sim->numberOfPointSets(); pointSetIndex_rk++) { + if (pointSetIndex_rk != pointSetIndex_r) { + Vector3r v_b_k_body = Vector3r::Zero(); + Vector3r omega_b_k_body = Vector3r::Zero(); + BoundaryModel_Akinci2012* bm_rk = static_cast(sim->getBoundaryModelFromPointSet(pointSetIndex_rk)); + int index_rk = pointSetIndex_rk - nFluids; + DynamicRigidBody* rb_rk = static_cast(bm_rk->getRigidBodyObject()); + for (unsigned int n = 0; n < sim->numberOfNeighbors(pointSetIndex_r, pointSetIndex_rk, r); n++) { + const unsigned int rk = sim->getNeighbor(pointSetIndex_r, pointSetIndex_rk, r, n); + Vector3r pos_rk = bm_rk->getPosition(rk) - rb_rk->getPosition(); + Vector3r grad_p_b_rkr = getDensity(index_rk, rk) * getArtificialVolume(index_r, r) * getDensity(index_r, r) / density_r2 * gradW(bm_rk->getPosition(rk) - bm_r->getPosition(r)); + v_b_k_body += -dt * rb_rk->getInvMass() * getArtificialVolume(index_rk, rk) * grad_p_b_rkr; + omega_b_k_body += -dt * rb_rk->getInertiaTensorInverseW() * getArtificialVolume(index_rk, rk) * pos_rk.cross(grad_p_b_rkr); + } + // divergence + Real sum_rk = 0; + for (unsigned int n = 0; n < sim->numberOfNeighbors(pointSetIndex_r, pointSetIndex_rk, r); n++) { + const unsigned int rk = sim->getNeighbor(pointSetIndex_r, pointSetIndex_rk, r, n); + Vector3r pos_rk = bm_rk->getPosition(rk) - rb_rk->getPosition(); + Vector3r v_b_rk = v_b_k_body + omega_b_k_body.cross(pos_rk); + sum_rk += getArtificialVolume(index_rk, rk) * getDensity(index_rk, rk) * (v_b_rk - v_b_r).dot(gradW(bm_r->getPosition(r) - bm_rk->getPosition(rk))); + } + b_r += sum_rk; + } + } + b_r = -b_r; + //setDiagonalElement(index_r, r, b_r < 0 ? b_r : 0); + setDiagonalElement(index_r, r, b_r); + } + } + } + } + } +} + +void StrongCouplingBoundarySolver::computeSourceTermRHS() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + computeSourceTermRHSForBody(boundaryPointSetIndex); + } +} + +void SPH::StrongCouplingBoundarySolver::computeSourceTermRHSForBody(const unsigned int& boundaryPointSetIndex) { + Simulation* sim = Simulation::getCurrent(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + + // compute pressure gradient, v_rr and omega_rr for rigid bodies + int bmIndex = boundaryPointSetIndex - nFluids; + DynamicRigidBody* rb = dynamic_cast(bm->getRigidBodyObject()); + Real vx = 0; + Real vy = 0; + Real vz = 0; + Real omegax = 0; + Real omegay = 0; + Real omegaz = 0; + if (rb != nullptr && rb->isDynamic() && m_bodyContacts[bmIndex] > 0) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) reduction(+:vx) reduction(+:vy) reduction(+:vz) reduction(+:omegax) reduction(+:omegay) reduction(+:omegaz) + for (int r = 0; r < bm->numberOfParticles(); r++) { + Vector3r pressureGrad_r = Vector3r::Zero(); + const Real density_r = getDensity(bmIndex, r); + const Real pressure_r = getPressure(bmIndex, r); + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + if (pid != boundaryPointSetIndex) { + BoundaryModel_Akinci2012* bm_neighbor = static_cast(sim->getBoundaryModelFromPointSet(pid)); + int neighborIndex = pid - nFluids; + for (unsigned int j = 0; j < sim->numberOfNeighbors(boundaryPointSetIndex, pid, r); j++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, pid, r, j); + const Real density_k = getDensity(neighborIndex, k); + const Real volume_k = getArtificialVolume(neighborIndex, k); + const Real pressure_k = getPressure(neighborIndex, k); + pressureGrad_r += volume_k * density_k * (pressure_r / (density_r * density_r) + pressure_k / (density_k * density_k)) * gradW(bm->getPosition(r) - bm_neighbor->getPosition(k)); + } + } + } + pressureGrad_r *= density_r; + setPressureGrad(bmIndex, r, pressureGrad_r); + + Vector3r v_rr_body_tmp = -dt * rb->getInvMass() * getArtificialVolume(bmIndex, r) * getPressureGrad(bmIndex, r); + vx += v_rr_body_tmp.x(); + vy += v_rr_body_tmp.y(); + vz += v_rr_body_tmp.z(); + Vector3r pos = bm->getPosition(r) - rb->getPosition(); + Vector3r omega_rr_body_tmp = -dt * rb->getInertiaTensorInverseW() * getArtificialVolume(bmIndex, r) * pos.cross(getPressureGrad(bmIndex, r)); + omegax += omega_rr_body_tmp.x(); + omegay += omega_rr_body_tmp.y(); + omegaz += omega_rr_body_tmp.z(); + } + } + } + Vector3r v_rr_body = Vector3r(vx, vy, vz); + Vector3r omega_rr_body = Vector3r(omegax, omegay, omegaz); + setV_rr_body(bmIndex, v_rr_body); + setOmega_rr_body(bmIndex, omega_rr_body); + + // compute v_rr and predicted velocity (v_rr+v_s) for all particles + if (rb != nullptr && rb->isDynamic() && m_bodyContacts[bmIndex] > 0) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm->numberOfParticles(); r++) { + Vector3r pos = bm->getPosition(r) - rb->getPosition(); + setV_rr(bmIndex, r, getV_rr_body(bmIndex) + getOmega_rr_body(bmIndex).cross(pos)); + setPredictedVelocity(bmIndex, r, getV_rr(bmIndex, r) + getV_s(bmIndex, r)); + setPredictedPosition(bmIndex, r, bm->getPosition(r) + (bm->getVelocity(r) + getPredictedVelocity(bmIndex, r)) / 2 * dt); + + Real minus_rho_div_v_rr = 0.0; + const Vector3r v_rr_r = getV_rr(bmIndex, r); + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + // divergence of particles in the same body should be 0; + if (boundaryPointSetIndex != pid) { + BoundaryModel_Akinci2012* bm_neighbor = static_cast(sim->getBoundaryModelFromPointSet(pid)); + int neighborIndex = pid - nFluids; + for (unsigned int j = 0; j < sim->numberOfNeighbors(boundaryPointSetIndex, pid, r); j++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, pid, r, j); + const Real density_k = getDensity(neighborIndex, k); + const Real volume_k = getArtificialVolume(neighborIndex, k); + const Vector3r v_rr_k = getV_rr(pid - nFluids, k); + minus_rho_div_v_rr += volume_k * density_k * (v_rr_k - v_rr_r).dot(gradW(bm->getPosition(r) - bm_neighbor->getPosition(k))); + } + } + } + minus_rho_div_v_rr = -minus_rho_div_v_rr; + setSourceTermRHS(bmIndex, r, minus_rho_div_v_rr); + } + } + } +} + + +void StrongCouplingBoundarySolver::pressureSolveIteration(Real& avgDensityDeviation) { + Simulation* sim = Simulation::getCurrent(); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + Real avgDensityDev = 0.0; + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + int bmIndex = boundaryPointSetIndex - nFluids; + if (bm->getRigidBodyObject()->isDynamic() && m_bodyContacts[bmIndex] > 0) { + int numContactsBody = getBodyContacts(bmIndex); + computeSourceTermRHSForBody(boundaryPointSetIndex); + // beta_r_RJ + Real relaxation = 0.5 / numContactsBody; + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) reduction(+:avgDensityDev) + for (int r = 0; r < bm->numberOfParticles(); r++) { + if (std::abs(getDiagonalElement(bmIndex, r)) > 1e-10) { + Real pressureNextIter = getPressure(bmIndex, r) + relaxation / getDiagonalElement(bmIndex, r) * (getSourceTerm(bmIndex, r) - getSourceTermRHS(bmIndex, r)); + setPressure(bmIndex, r, std::max(pressureNextIter, static_cast(0.0))); + avgDensityDev -= (getSourceTerm(bmIndex, r) - getMinus_rho_div_v_s(bmIndex, r) - getSourceTermRHS(bmIndex, r)) * dt; + } else { + setPressure(bmIndex, r, static_cast(0.0)); + } + } + } + } + } + avgDensityDev /= m_contactsAllBodies; + avgDensityDev /= getRestDensity(); + avgDensityDeviation = std::abs(avgDensityDev); +} + +void StrongCouplingBoundarySolver::applyForce() { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + for (unsigned int boundaryPointSetIndex = nFluids; boundaryPointSetIndex < sim->numberOfPointSets(); boundaryPointSetIndex++) { + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + int bmIndex = boundaryPointSetIndex - nFluids; + DynamicRigidBody* rb = dynamic_cast(bm->getRigidBodyObject()); + if (rb != nullptr && rb->isDynamic()) { + #pragma omp parallel default(shared) + { + #pragma omp for schedule(static) + for (int r = 0; r < bm->numberOfParticles(); r++) { + const Vector3r f = -getArtificialVolume(bmIndex, r) * getPressureGrad(bmIndex, r); + bm->addForce(bm->getPosition(r), f); + if (f.norm() > 1e-10) { + const Vector3r fric = copmuteParticleFriction(boundaryPointSetIndex, r, f); + bm->addForce(bm->getPosition(r), fric); + } + setPressure(bmIndex, r, static_cast(0.0)); + setPressureGrad(bmIndex, r, Vector3r::Zero()); + } + } + } + } +} + +Vector3r StrongCouplingBoundarySolver::copmuteParticleFriction(const unsigned int& boundaryPointSetIndex, const unsigned int& index, const Vector3r& force) { + Simulation* sim = Simulation::getCurrent(); + DynamicBoundarySimulator* boundarySimulator = sim->getDynamicBoundarySimulator(); + TimeManager* tm = TimeManager::getCurrent(); + const Real dt = tm->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + BoundaryModel_Akinci2012* bm = static_cast(sim->getBoundaryModelFromPointSet(boundaryPointSetIndex)); + DynamicRigidBody* rb = dynamic_cast(bm->getRigidBodyObject()); + if (rb == nullptr) { + return Vector3r::Zero(); + } + int bmIndex = boundaryPointSetIndex - nFluids; + + // compute normal + Vector3r numerator = bm->getPosition(index) * W_zero(); + Real denom = W_zero(); + for (unsigned int n = 0; n < sim->numberOfNeighbors(boundaryPointSetIndex, boundaryPointSetIndex, index); n++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, boundaryPointSetIndex, index, n); + numerator += bm->getPosition(k) * W(bm->getPosition(index) - bm->getPosition(k)); + denom += W(bm->getPosition(index) - bm->getPosition(k)); + } + Vector3r normal = bm->getPosition(index) - numerator / denom; + normal.normalize(); + + // compute averaged relative predicted velocity + Vector3r numeratorV = Vector3r::Zero(); + Real denomV = 0; + for (unsigned int pid = nFluids; pid < sim->numberOfPointSets(); pid++) { + BoundaryModel_Akinci2012* bm_k = static_cast(sim->getBoundaryModelFromPointSet(pid)); + if (pid != boundaryPointSetIndex) { + for (unsigned int n = 0; n < sim->numberOfNeighbors(boundaryPointSetIndex, pid, index); n++) { + const unsigned int k = sim->getNeighbor(boundaryPointSetIndex, pid, index, n); + numeratorV += getPredictedVelocity(pid - nFluids, k) * W(bm->getPosition(index) - bm_k->getPosition(k)); + denomV += W(bm->getPosition(index) - bm_k->getPosition(k)); + } + } + } + Vector3r neighborV = numeratorV / denomV; + Vector3r relativeV = getPredictedVelocity(bmIndex, index) - neighborV; + Vector3r normalVPred = relativeV.dot(normal) * normal; + Vector3r tangentialVPred = relativeV - normalVPred; + + Vector3r tangent = tangentialVPred; + tangent.normalize(); + Vector3r normalF = -force.dot(normal) * normal; + Vector3r tangentF = -rb->getFrictionCoeff() * normalF.norm() * tangent; + Vector3r pos = bm->getPosition(index); + Matrix3r productMat; + productMat << 0, -pos.z(), pos.y(), + pos.z(), 0, -pos.x(), + -pos.y(), pos.x(), 0; + Matrix3r K = rb->getInvMass() * Matrix3r::Identity() - productMat * rb->getInertiaTensorInverseW() * productMat; + Vector3r max = -tangentialVPred / (dt * getBodyContacts(bmIndex) * tangent.transpose() * K * tangent); + if (max.norm() < tangentF.norm()) { + tangentF = max; + } + return tangentF; +} + +Real StrongCouplingBoundarySolver::W(const Vector3r& r) { + return m_kernelFct(r); +} + +Vector3r StrongCouplingBoundarySolver::gradW(const Vector3r& r) { + return m_gradKernelFct(r); +} + +Real StrongCouplingBoundarySolver::W_zero() { + return m_W_zero; +} + +void StrongCouplingBoundarySolver::setKernel(int val) { + Simulation* sim = Simulation::getCurrent(); + if (val == m_kernelMethod) + return; + + m_kernelMethod = val; + if (!sim->is2DSimulation()) { + if ((m_kernelMethod < 0) || (m_kernelMethod > 4)) + m_kernelMethod = 0; + + if (m_kernelMethod == 0) { + m_W_zero = CubicKernel::W_zero(); + m_kernelFct = CubicKernel::W; + } else if (m_kernelMethod == 1) { + m_W_zero = WendlandQuinticC2Kernel::W_zero(); + m_kernelFct = WendlandQuinticC2Kernel::W; + } else if (m_kernelMethod == 2) { + m_W_zero = Poly6Kernel::W_zero(); + m_kernelFct = Poly6Kernel::W; + } else if (m_kernelMethod == 3) { + m_W_zero = SpikyKernel::W_zero(); + m_kernelFct = SpikyKernel::W; + } else if (m_kernelMethod == 4) { + m_W_zero = Simulation::PrecomputedCubicKernel::W_zero(); + m_kernelFct = Simulation::PrecomputedCubicKernel::W; + } + } else { + if ((m_kernelMethod < 0) || (m_kernelMethod > 1)) + m_kernelMethod = 0; + + if (m_kernelMethod == 0) { + m_W_zero = CubicKernel2D::W_zero(); + m_kernelFct = CubicKernel2D::W; + } else if (m_kernelMethod == 1) { + m_W_zero = WendlandQuinticC2Kernel2D::W_zero(); + m_kernelFct = WendlandQuinticC2Kernel2D::W; + } + } + + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + sim->updateBoundaryVolume(); +} + +void StrongCouplingBoundarySolver::setGradKernel(int val) { + Simulation* sim = Simulation::getCurrent(); + m_gradKernelMethod = val; + if (!sim->is2DSimulation()) { + if ((m_gradKernelMethod < 0) || (m_gradKernelMethod > 4)) + m_gradKernelMethod = 0; + + if (m_gradKernelMethod == 0) + m_gradKernelFct = CubicKernel::gradW; + else if (m_gradKernelMethod == 1) + m_gradKernelFct = WendlandQuinticC2Kernel::gradW; + else if (m_gradKernelMethod == 2) + m_gradKernelFct = Poly6Kernel::gradW; + else if (m_gradKernelMethod == 3) + m_gradKernelFct = SpikyKernel::gradW; + else if (m_gradKernelMethod == 4) + m_gradKernelFct = Simulation::PrecomputedCubicKernel::gradW; + } else { + if ((m_gradKernelMethod < 0) || (m_gradKernelMethod > 1)) + m_gradKernelMethod = 0; + + if (m_gradKernelMethod == 0) + m_gradKernelFct = CubicKernel2D::gradW; + else if (m_gradKernelMethod == 1) + m_gradKernelFct = WendlandQuinticC2Kernel2D::gradW; + } +} diff --git a/SPlisHSPlasH/StrongCouplingBoundarySolver.h b/SPlisHSPlasH/StrongCouplingBoundarySolver.h new file mode 100644 index 00000000..717e8dd9 --- /dev/null +++ b/SPlisHSPlasH/StrongCouplingBoundarySolver.h @@ -0,0 +1,321 @@ +#include "SPlisHSPlasH/Common.h" +#include "SPlisHSPlasH/TimeStep.h" +#include "SPlisHSPlasH/SPHKernels.h" +#include "ParameterObject.h" +#include + +namespace SPH { + /** \brief This class stores the information of a dynamic rigid body which + * is used for the strong coupling method introduced in + * Interlinked SPH Pressure Solvers for Strong Fluid-Rigid Coupling. Gissler et al. https://doi.org/10.1145/3284980 + */ + class StrongCouplingBoundarySolver : public GenParam::ParameterObject { + private: + // values required for Gissler 2019 strong coupling based on Akinci 2012 + Real m_restDensity; + std::vector> m_density; + std::vector> m_pressure; + std::vector> m_artificialVolume; + std::vector> m_v_s; + std::vector> m_s; // source term + std::vector> m_pressureGrad; + std::vector> m_v_rr; + std::vector> m_predictedVelocity; + std::vector> m_predictedPosition; + std::vector> m_minus_rho_div_v_s; + std::vector> m_minus_rho_div_v_rr; // RHS to the source term + std::vector> m_diagonalElement; // diagonal element for jacobi iteration + std::vector m_v_rr_body; + std::vector m_omega_rr_body; + std::vector m_bodyContacts; + unsigned int m_contactsAllBodies; + + static StrongCouplingBoundarySolver* current; + + // only used in WCSPH + unsigned int m_maxIterations; + unsigned int m_minIterations; + + Real m_maxDensityDeviation; + + int m_kernelMethod; + int m_gradKernelMethod; + Real(*m_kernelFct)(const Vector3r&); + Vector3r(*m_gradKernelFct)(const Vector3r& r); + Real m_W_zero; + + public: + static int KERNEL_METHOD; + static int GRAD_KERNEL_METHOD; + static int ENUM_KERNEL_CUBIC; + static int ENUM_KERNEL_WENDLANDQUINTICC2; + static int ENUM_KERNEL_POLY6; + static int ENUM_KERNEL_SPIKY; + static int ENUM_KERNEL_CUBIC_2D; + static int ENUM_KERNEL_WENDLANDQUINTICC2_2D; + static int ENUM_GRADKERNEL_CUBIC; + static int ENUM_GRADKERNEL_WENDLANDQUINTICC2; + static int ENUM_GRADKERNEL_POLY6; + static int ENUM_GRADKERNEL_SPIKY; + static int ENUM_GRADKERNEL_CUBIC_2D; + static int ENUM_GRADKERNEL_WENDLANDQUINTICC2_2D; + + StrongCouplingBoundarySolver(); + ~StrongCouplingBoundarySolver(); + static StrongCouplingBoundarySolver* getCurrent(); + virtual void initParameters(); + void reset(); + void resize(unsigned int size); + void computeContacts(); + void computeDensityAndVolume(); + void computeV_s(); + void computeSourceTermRHS(); + void computeSourceTermRHSForBody(const unsigned int& boundaryPointSetIndex); + void computeSourceTerm(); + void computeDiagonalElement(); + void pressureSolveIteration(Real& avgDensityDeviation); + Vector3r copmuteParticleFriction(const unsigned int& boundaryPointSetIndex, const unsigned int& index, const Vector3r& force); + void applyForce(); + void performNeighborhoodSearchSort(); + Real W(const Vector3r& r); + Vector3r gradW(const Vector3r& r); + Real W_zero(); + + + + FORCE_INLINE const Real& getDensity(const unsigned int& rigidBodyIndex, const unsigned int& index) const { + return m_density[rigidBodyIndex][index]; + } + + FORCE_INLINE Real& getDensity(const unsigned int& rigidBodyIndex, const unsigned int& index) { + return m_density[rigidBodyIndex][index]; + } + + FORCE_INLINE void setDensity(const unsigned int& rigidBodyIndex, const unsigned int& index, const Real& value) { + m_density[rigidBodyIndex][index] = value; + } + + FORCE_INLINE const Real& getPressure(const unsigned int& rigidBodyIndex, const unsigned int& index) const { + return m_pressure[rigidBodyIndex][index]; + } + + FORCE_INLINE Real& getPressure(const unsigned int& rigidBodyIndex, const unsigned int& index) { + return m_pressure[rigidBodyIndex][index]; + } + + FORCE_INLINE void setPressure(const unsigned int& rigidBodyIndex, const unsigned int& index, const Real& value) { + m_pressure[rigidBodyIndex][index] = value; + } + + FORCE_INLINE const Real& getRestDensity() const { + return m_restDensity; + } + + FORCE_INLINE Real& getRestDensity() { + return m_restDensity; + } + + FORCE_INLINE void setRestDensity(const Real& value) { + m_restDensity = value; + } + + FORCE_INLINE const Vector3r& getV_rr_body(const unsigned int& rigidBodyIndex) const { + return m_v_rr_body[rigidBodyIndex]; + } + + FORCE_INLINE Vector3r& getV_rr_body(const unsigned int& rigidBodyIndex) { + return m_v_rr_body[rigidBodyIndex]; + } + + FORCE_INLINE void setV_rr_body(const unsigned int& rigidBodyIndex, const Vector3r& value) { + m_v_rr_body[rigidBodyIndex] = value; + } + + FORCE_INLINE const Vector3r& getOmega_rr_body(const unsigned int& rigidBodyIndex) const { + return m_omega_rr_body[rigidBodyIndex]; + } + + FORCE_INLINE Vector3r& getOmega_rr_body(const unsigned int& rigidBodyIndex) { + return m_omega_rr_body[rigidBodyIndex]; + } + + FORCE_INLINE void setOmega_rr_body(const unsigned int& rigidBodyIndex, const Vector3r& value) { + m_omega_rr_body[rigidBodyIndex] = value; + } + + FORCE_INLINE Vector3r& getV_s(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_v_s[rigidBodyIndex][i]; + } + + FORCE_INLINE const Vector3r& getV_s(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_v_s[rigidBodyIndex][i]; + } + + FORCE_INLINE void setV_s(const unsigned int& rigidBodyIndex, const unsigned int i, const Vector3r& value) { + m_v_s[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Vector3r& getV_rr(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_v_rr[rigidBodyIndex][i]; + } + + FORCE_INLINE const Vector3r& getV_rr(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_v_rr[rigidBodyIndex][i]; + } + + FORCE_INLINE void setV_rr(const unsigned int& rigidBodyIndex, const unsigned int i, const Vector3r& value) { + m_v_rr[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Vector3r& getPredictedVelocity(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_predictedVelocity[rigidBodyIndex][i]; + } + + FORCE_INLINE const Vector3r& getPredictedVelocity(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_predictedVelocity[rigidBodyIndex][i]; + } + + FORCE_INLINE void setPredictedVelocity(const unsigned int& rigidBodyIndex, const unsigned int i, const Vector3r& value) { + m_predictedVelocity[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Vector3r& getPredictedPosition(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_predictedPosition[rigidBodyIndex][i]; + } + + FORCE_INLINE const Vector3r& getPredictedPosition(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_predictedPosition[rigidBodyIndex][i]; + } + + FORCE_INLINE void setPredictedPosition(const unsigned int& rigidBodyIndex, const unsigned int i, const Vector3r& value) { + m_predictedPosition[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Vector3r& getPressureGrad(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_pressureGrad[rigidBodyIndex][i]; + } + + FORCE_INLINE const Vector3r& getPressureGrad(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_pressureGrad[rigidBodyIndex][i]; + } + + FORCE_INLINE void setPressureGrad(const unsigned int& rigidBodyIndex, const unsigned int i, const Vector3r& value) { + m_pressureGrad[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Real& getMinus_rho_div_v_s(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_minus_rho_div_v_s[rigidBodyIndex][i]; + } + + FORCE_INLINE const Real& getMinus_rho_div_v_s(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_minus_rho_div_v_s[rigidBodyIndex][i]; + } + + FORCE_INLINE void setMinus_rho_div_v_s(const unsigned int& rigidBodyIndex, const unsigned int i, const Real& value) { + m_minus_rho_div_v_s[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Real& getSourceTermRHS(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_minus_rho_div_v_rr[rigidBodyIndex][i]; + } + + FORCE_INLINE const Real& getSourceTermRHS(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_minus_rho_div_v_rr[rigidBodyIndex][i]; + } + + FORCE_INLINE void setSourceTermRHS(const unsigned int& rigidBodyIndex, const unsigned int i, const Real& value) { + m_minus_rho_div_v_rr[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Real& getDiagonalElement(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_diagonalElement[rigidBodyIndex][i]; + } + + FORCE_INLINE const Real& getDiagonalElement(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_diagonalElement[rigidBodyIndex][i]; + } + + FORCE_INLINE void setDiagonalElement(const unsigned int& rigidBodyIndex, const unsigned int i, const Real& value) { + m_diagonalElement[rigidBodyIndex][i] = value; + } + + FORCE_INLINE Real& getSourceTerm(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_s[rigidBodyIndex][i]; + } + + FORCE_INLINE const Real& getSourceTerm(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_s[rigidBodyIndex][i]; + } + + FORCE_INLINE void setSourceTerm(const unsigned int& rigidBodyIndex, const unsigned int i, const Real& value) { + m_s[rigidBodyIndex][i] = value; + } + + FORCE_INLINE const Real& getArtificialVolume(const unsigned int& rigidBodyIndex, const unsigned int i) const { + return m_artificialVolume[rigidBodyIndex][i]; + } + + FORCE_INLINE Real& getArtificialVolume(const unsigned int& rigidBodyIndex, const unsigned int i) { + return m_artificialVolume[rigidBodyIndex][i]; + } + + FORCE_INLINE void setArtificialVolume(const unsigned int& rigidBodyIndex, const unsigned int i, const Real& val) { + m_artificialVolume[rigidBodyIndex][i] = val; + } + + FORCE_INLINE const unsigned int& getBodyContacts(const unsigned int& rigidBodyIndex) const { + return m_bodyContacts[rigidBodyIndex]; + } + + FORCE_INLINE unsigned int& getBodyContacts(const unsigned int& rigidBodyIndex) { + return m_bodyContacts[rigidBodyIndex]; + } + + FORCE_INLINE void setBodyContacts(const unsigned int& rigidBodyIndex, const unsigned int& value) { + m_bodyContacts[rigidBodyIndex] = value; + } + + FORCE_INLINE unsigned int getAllContacts() { + return m_contactsAllBodies; + } + + FORCE_INLINE const unsigned int& getMaxIterations() const { + return m_maxIterations; + } + FORCE_INLINE unsigned int& getMaxIterations() { + return m_maxIterations; + } + FORCE_INLINE void setMaxIterations(const unsigned int& value) { + m_maxIterations = value; + } + + FORCE_INLINE const unsigned int& getMinIterations() const { + return m_minIterations; + } + FORCE_INLINE unsigned int& getMinIterations() { + return m_minIterations; + } + FORCE_INLINE void setMinIterations(const unsigned int& value) { + m_minIterations = value; + } + + FORCE_INLINE const Real& getMaxDensityDeviation() const { + return m_maxDensityDeviation; + } + FORCE_INLINE Real& getMaxDensityDeviation() { + return m_maxDensityDeviation; + } + FORCE_INLINE void setMaxDensityDeviation(const Real& value) { + m_maxDensityDeviation = value; + } + + int getKernel() const { + return m_kernelMethod; + } + void setKernel(int val); + int getGradKernel() const { + return m_gradKernelMethod; + } + void setGradKernel(int val); + }; +} \ No newline at end of file diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp index 867919c5..d531679a 100644 --- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp +++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp @@ -109,7 +109,7 @@ void SurfaceTension_Akinci2013::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) { forall_boundary_neighbors( // adhesion force diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp index 303d556d..326d4742 100644 --- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp +++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp @@ -56,7 +56,7 @@ void SurfaceTension_Becker2007::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r xixj = xi - xj; diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp index e9187385..01142b24 100644 --- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp +++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp @@ -63,7 +63,7 @@ void SurfaceTension_He2014::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( ci += bm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj); @@ -132,7 +132,7 @@ void SurfaceTension_He2014::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r a = factorb*bm_neighbor->getVolume(neighborIndex) * gradC2_i * sim->gradW(xi - xj); diff --git a/SPlisHSPlasH/TimeStep.cpp b/SPlisHSPlasH/TimeStep.cpp index e4353878..b65e27ed 100644 --- a/SPlisHSPlasH/TimeStep.cpp +++ b/SPlisHSPlasH/TimeStep.cpp @@ -113,7 +113,7 @@ void TimeStep::computeDensities(const unsigned int fluidModelIndex) ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Scalarf8 V_avx = convert_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVolume(0), count); diff --git a/SPlisHSPlasH/Utilities/SceneLoader.cpp b/SPlisHSPlasH/Utilities/SceneLoader.cpp index d1c55150..b9b31bb0 100644 --- a/SPlisHSPlasH/Utilities/SceneLoader.cpp +++ b/SPlisHSPlasH/Utilities/SceneLoader.cpp @@ -44,6 +44,9 @@ void SceneLoader::readScene(const char *fileName, Scene &scene) scene.sim2D = false; readValue(config["sim2D"], scene.sim2D); + + scene.boundaryHandlingMethod = 0; + readValue(config["boundaryHandlingMethod"], scene.boundaryHandlingMethod); } ////////////////////////////////////////////////////////////////////////// diff --git a/SPlisHSPlasH/Utilities/SceneLoader.h b/SPlisHSPlasH/Utilities/SceneLoader.h index a404bd9e..9cceb17c 100644 --- a/SPlisHSPlasH/Utilities/SceneLoader.h +++ b/SPlisHSPlasH/Utilities/SceneLoader.h @@ -30,6 +30,7 @@ namespace Utilities std::vector materials; Real particleRadius; bool sim2D; + int boundaryHandlingMethod; }; nlohmann::json& getJSONData() { return m_jsonData; }; diff --git a/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp b/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp index 87431e97..ac05a172 100644 --- a/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp +++ b/SPlisHSPlasH/Utilities/SceneParameterObjects.cpp @@ -318,6 +318,9 @@ int BoundaryParameterObject::BOUNDARY_MAP_THICKNESS = -1; int BoundaryParameterObject::BOUNDARY_MAP_RESOLUTION = -1; int BoundaryParameterObject::BOUNDARY_SAMPLING_MODE = -1; int BoundaryParameterObject::BOUNDARY_IS_ANIMATED = -1; +int BoundaryParameterObject::BOUNDARY_DENSITY = -1; +int BoundaryParameterObject::BOUNDARY_VELOCITY = -1; +int BoundaryParameterObject::BOUNDARY_FRICTION = -1; void BoundaryParameterObject::initParameters() { @@ -380,4 +383,16 @@ void BoundaryParameterObject::initParameters() BOUNDARY_IS_ANIMATED = createBoolParameter("isAnimated", "Animated", &isAnimated); setGroup(BOUNDARY_IS_ANIMATED, "Boundary"); setDescription(BOUNDARY_IS_ANIMATED, "Defines if the body is animated (e.g. by a script)."); + + BOUNDARY_DENSITY = createNumericParameter("density", "Density", &density); + setGroup(BOUNDARY_DENSITY, "Boundary"); + setDescription(BOUNDARY_DENSITY, "Density of the rigid body."); + + BOUNDARY_FRICTION = createNumericParameter("friction", "Friction", &friction); + setGroup(BOUNDARY_FRICTION, "Boundary"); + setDescription(BOUNDARY_FRICTION, "Friction of the rigid body."); + + BOUNDARY_VELOCITY = createVectorParameter("velocity", "Velocity", 3u, velocity.data()); + setGroup(BOUNDARY_VELOCITY, "Boundary"); + setDescription(BOUNDARY_VELOCITY, "Velocity of the rigid body after loading."); } \ No newline at end of file diff --git a/SPlisHSPlasH/Utilities/SceneParameterObjects.h b/SPlisHSPlasH/Utilities/SceneParameterObjects.h index e2bba70e..b9957c53 100644 --- a/SPlisHSPlasH/Utilities/SceneParameterObjects.h +++ b/SPlisHSPlasH/Utilities/SceneParameterObjects.h @@ -322,8 +322,11 @@ namespace Utilities std::string meshFile; Vector3r translation; Vector3r axis; + Vector3r velocity; Real angle; Vector3r scale; + Real density; + Real friction; bool dynamic; bool isWall; Eigen::Matrix color; @@ -341,10 +344,13 @@ namespace Utilities meshFile = ""; translation = Vector3r::Zero(); axis = Vector3r(1, 0, 0); + velocity = Vector3r::Zero(); angle = 0.0; + density = 1000.0; scale = Vector3r::Ones(); dynamic = false; isWall = false; + friction = 0.2; color = Eigen::Vector4f(1.0f, 0.0f, 0.0f, 0.0f); samplingMode = 0; isAnimated = false; @@ -357,15 +363,19 @@ namespace Utilities BoundaryParameterObject(std::string samplesFile_, std::string meshFile_, Vector3r translation_, Vector3r axis_, Real angle_, Vector3r scale_, bool dynamic_, bool isWall_, Eigen::Matrix color_, std::string mapFile_, bool mapInvert_, - Real mapThickness_, Eigen::Matrix mapResolution_, unsigned int samplingMode_, bool isAnimated_) + Real mapThickness_, Eigen::Matrix mapResolution_, unsigned int samplingMode_, bool isAnimated_, + Vector3r velocity_, Real friction_, Real density_ = 1000.0) { samplesFile = samplesFile_; meshFile = meshFile_; translation = translation_; axis = axis_; angle = angle_; + velocity = velocity_; scale = scale_; + density = density_; dynamic = dynamic_; + friction = friction_; isWall = isWall_; color = color_; samplingMode = samplingMode_; @@ -392,6 +402,9 @@ namespace Utilities static int BOUNDARY_MAP_RESOLUTION; static int BOUNDARY_SAMPLING_MODE; static int BOUNDARY_IS_ANIMATED; + static int BOUNDARY_DENSITY; + static int BOUNDARY_VELOCITY; + static int BOUNDARY_FRICTION; virtual void initParameters(); }; diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp index 22ad1954..a78c2fb7 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp @@ -201,7 +201,7 @@ void Viscosity_Bender2017::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp index fb51e952..ca098070 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp @@ -89,7 +89,7 @@ void Viscosity_Peer2015::computeDensities() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( // Boundary: Akinci2012 diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp index 64665ac7..96647c58 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp @@ -114,7 +114,7 @@ void Viscosity_Peer2016::computeDensities() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( // Boundary: Akinci2012 diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp index f91671ff..1fec59e1 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp @@ -94,7 +94,7 @@ void Viscosity_Standard::step() ////////////////////////////////////////////////////////////////////////// if (m_boundaryViscosity != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVelocity(0), count); diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp index 9a72f5e5..3b9425e8 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp @@ -214,7 +214,7 @@ void Viscosity_Takahashi2015::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp index 56469502..83279e03 100644 --- a/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp +++ b/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp @@ -130,7 +130,7 @@ void Viscosity_Weiler2018::matrixVecProd(const Real* vec, Real *result, void *us ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVelocity(0), count); @@ -492,7 +492,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Matrix3r ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Vector3f8 xixj = xi_avx - xj_avx; @@ -636,7 +636,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Matrix3r ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r xixj = xi - xj; @@ -779,7 +779,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Vector3r ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r xixj = xi - xj; @@ -954,7 +954,7 @@ void Viscosity_Weiler2018::applyForces(const VectorXr &x) ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); @@ -1118,7 +1118,7 @@ void Viscosity_Weiler2018::computeRHS(VectorXr &b, VectorXr &g) ////////////////////////////////////////////////////////////////////////// if (mub != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); diff --git a/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp b/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp index 52f60465..c6919f25 100644 --- a/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp +++ b/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp @@ -135,7 +135,7 @@ void MicropolarModel_Bender2017::step() ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors_avx( const Vector3f8 vj_avx = convertVec_zero(&sim->getNeighborList(fluidModelIndex, pid, i)[j], &bm_neighbor->getVelocity(0), count); diff --git a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp index 6fd56500..4abb4236 100644 --- a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp +++ b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp @@ -8,7 +8,9 @@ #include "SPlisHSPlasH/BoundaryModel_Akinci2012.h" #include "SPlisHSPlasH/BoundaryModel_Koschier2017.h" #include "SPlisHSPlasH/BoundaryModel_Bender2019.h" - +#include "Simulator/DynamicBoundarySimulator.h" +#include "SPlisHSPlasH/DynamicRigidBody.h" +#include "SPlisHSPlasH/StrongCouplingBoundarySolver.h" using namespace SPH; using namespace std; @@ -91,6 +93,7 @@ void TimeStepWCSPH::step() } sim->computeNonPressureForces(); + sim->updateTimeStepSize(); for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++) { @@ -110,7 +113,9 @@ void TimeStepWCSPH::step() computePressureAccels(fluidModelIndex); } - sim->updateTimeStepSize(); + if (sim->getDynamicBoundarySimulator() != nullptr && sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + solveRigidRigidContacts(); + } for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++) { @@ -137,81 +142,99 @@ void TimeStepWCSPH::step() sim->animateParticles(); // Compute new time - tm->setTime (tm->getTime () + h); + tm->setTime(tm->getTime() + h); } -void TimeStepWCSPH::reset() -{ +void TimeStepWCSPH::reset() { TimeStep::reset(); m_simulationData.reset(); m_counter = 0; } -void TimeStepWCSPH::computePressureAccels(const unsigned int fluidModelIndex) -{ - Simulation *sim = Simulation::getCurrent(); - FluidModel *model = sim->getFluidModel(fluidModelIndex); +void TimeStepWCSPH::computePressureAccels(const unsigned int fluidModelIndex) { + Simulation* sim = Simulation::getCurrent(); + FluidModel* model = sim->getFluidModel(fluidModelIndex); const Real density0 = model->getDensity0(); const unsigned int numParticles = model->numActiveParticles(); const unsigned int nFluids = sim->numberOfFluidModels(); const unsigned int nBoundaries = sim->numberOfBoundaryModels(); // Compute pressure forces - #pragma omp parallel default(shared) + #pragma omp parallel default(shared) { - #pragma omp for schedule(static) - for (int i = 0; i < (int)numParticles; i++) - { - const Vector3r &xi = model->getPosition(i); + #pragma omp for schedule(static) + for (int i = 0; i < (int)numParticles; i++) { + const Vector3r& xi = model->getPosition(i); const Real density_i = model->getDensity(i); - Vector3r &ai = m_simulationData.getPressureAccel(fluidModelIndex, i); + Vector3r& ai = m_simulationData.getPressureAccel(fluidModelIndex, i); ai.setZero(); - const Real dpi = m_simulationData.getPressure(fluidModelIndex, i) / (density_i*density_i); + const Real dpi = m_simulationData.getPressure(fluidModelIndex, i) / (density_i * density_i); ////////////////////////////////////////////////////////////////////////// - // Fluid + // Fluid ////////////////////////////////////////////////////////////////////////// forall_fluid_neighbors( // Pressure const Real density_j = fm_neighbor->getDensity(neighborIndex) * density0 / fm_neighbor->getDensity0(); - const Real dpj = m_simulationData.getPressure(pid, neighborIndex) / (density_j*density_j); - ai -= density0 * fm_neighbor->getVolume(neighborIndex) * (dpi + dpj) * sim->gradW(xi - xj); + const Real dpj = m_simulationData.getPressure(pid, neighborIndex) / (density_j * density_j); + ai -= density0 * fm_neighbor->getVolume(neighborIndex) * (dpi + dpj) * sim->gradW(xi - xj); ); ////////////////////////////////////////////////////////////////////////// // Boundary ////////////////////////////////////////////////////////////////////////// - const Real dpj = m_simulationData.getPressure(fluidModelIndex, i) / (density0*density0); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) - { + const Real dpj = m_simulationData.getPressure(fluidModelIndex, i) / (density0 * density0); + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( - const Vector3r a = density0 * bm_neighbor->getVolume(neighborIndex) * (dpi + dpj)* sim->gradW(xi - xj); - ai -= a; - bm_neighbor->addForce(xj, model->getMass(i) * a); + const Vector3r a = density0 * bm_neighbor->getVolume(neighborIndex) * (dpi + dpj) * sim->gradW(xi - xj); + ai -= a; + bm_neighbor->addForce(xj, model->getMass(i) * a); ); - } - else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) - { + } else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) { forall_density_maps( - const Vector3r a = -density0 * (dpi + dpj)* gradRho; - ai -= a; - bm_neighbor->addForce(xj, model->getMass(i) * a); + const Vector3r a = -density0 * (dpi + dpj) * gradRho; + ai -= a; + bm_neighbor->addForce(xj, model->getMass(i) * a); ); - } - else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) - { + } else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) { forall_volume_maps( - const Vector3r a = density0 * Vj * (dpi + dpj)* sim->gradW(xi - xj); - ai -= a; - bm_neighbor->addForce(xj, model->getMass(i) * a); + const Vector3r a = density0 * Vj * (dpi + dpj) * sim->gradW(xi - xj); + ai -= a; + bm_neighbor->addForce(xj, model->getMass(i) * a); ); } } } } +void TimeStepWCSPH::solveRigidRigidContacts() { + Simulation* sim = Simulation::getCurrent(); + const Real dt = TimeManager::getCurrent()->getTimeStepSize(); + const unsigned int nFluids = sim->numberOfFluidModels(); + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + sim->getDynamicBoundarySimulator()->updateBoundaryForces(); + // check whehter there is any collisions + bs->computeContacts(); + if (bs->getAllContacts() == 0) { + return; + } + + bs->computeDensityAndVolume(); + bs->computeDiagonalElement(); + bs->computeSourceTerm(); + Real avgDensityDeviation = 1.0; + int iterations = 0; + // while density deviation too large + while ((avgDensityDeviation > bs->getMaxDensityDeviation() && iterations < bs->getMaxIterations()) || iterations < bs->getMinIterations()) { + avgDensityDeviation = 0.0; + bs->pressureSolveIteration(avgDensityDeviation); + iterations++; + } + bs->applyForce(); +} + void TimeStepWCSPH::performNeighborhoodSearch() { if (Simulation::getCurrent()->zSortEnabled()) diff --git a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.h b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.h index 598a0e6e..d1cfdd81 100644 --- a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.h +++ b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.h @@ -35,6 +35,8 @@ namespace SPH virtual void emittedParticles(FluidModel *model, const unsigned int startIndex); virtual void initParameters(); + void solveRigidRigidContacts(); + public: static int STIFFNESS; static int EXPONENT; diff --git a/SPlisHSPlasH/XSPH.cpp b/SPlisHSPlasH/XSPH.cpp index 1559161d..98e6f253 100644 --- a/SPlisHSPlasH/XSPH.cpp +++ b/SPlisHSPlasH/XSPH.cpp @@ -98,7 +98,7 @@ void XSPH::step() ////////////////////////////////////////////////////////////////////////// if (m_boundaryCoefficient != 0.0) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { forall_boundary_neighbors( const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex); diff --git a/Simulator/CMakeLists.txt b/Simulator/CMakeLists.txt index 2ff4af93..f354de07 100644 --- a/Simulator/CMakeLists.txt +++ b/Simulator/CMakeLists.txt @@ -56,6 +56,7 @@ set(EXPORTER_SOURCE_FILES ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_OBJ.cpp ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_PLY.cpp ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_VTK.cpp + ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyParticleExporter_VTK.cpp ) set(EXPORTER_HEADER_FILES @@ -66,6 +67,7 @@ set(EXPORTER_HEADER_FILES ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_OBJ.h ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_PLY.h ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyExporter_VTK.h + ${PROJECT_SOURCE_DIR}/Simulator/Exporter/RigidBodyParticleExporter_VTK.h ) if (USE_EMBEDDED_PYTHON) @@ -87,6 +89,8 @@ add_library(SimulatorBase ${PROJECT_SOURCE_DIR}/Simulator/BoundarySimulator.h ${PROJECT_SOURCE_DIR}/Simulator/StaticBoundarySimulator.cpp ${PROJECT_SOURCE_DIR}/Simulator/StaticBoundarySimulator.h + ${PROJECT_SOURCE_DIR}/Simulator/DynamicBoundarySimulator.cpp + ${PROJECT_SOURCE_DIR}/Simulator/DynamicBoundarySimulator.h ${PROJECT_SOURCE_DIR}/Simulator/GUI/OpenGL/Simulator_OpenGL.cpp ${PROJECT_SOURCE_DIR}/Simulator/GUI/OpenGL/Simulator_OpenGL.h diff --git a/Simulator/DynamicBoundarySimulator.cpp b/Simulator/DynamicBoundarySimulator.cpp new file mode 100644 index 00000000..b15b79ab --- /dev/null +++ b/Simulator/DynamicBoundarySimulator.cpp @@ -0,0 +1,232 @@ +#include "DynamicBoundarySimulator.h" +#include "SimulatorBase.h" +#include "Utilities/FileSystem.h" +#include "SPlisHSPlasH/Simulation.h" +#include "Utilities/PartioReaderWriter.h" +#include "SPlisHSPlasH/DynamicRigidBody.h" +#include "Utilities/Logger.h" +#include "Utilities/Timing.h" +#include "SPlisHSPlasH/Utilities/SurfaceSampling.h" +#include "SPlisHSPlasH/TriangleMesh.h" +#include "Simulator/SceneConfiguration.h" +#include "SPlisHSPlasH/Utilities/MeshImport.h" +#include "SPlisHSPlasH/StrongCouplingBoundarySolver.h" + +using namespace std; +using namespace SPH; +using namespace Utilities; + + +DynamicBoundarySimulator::DynamicBoundarySimulator(SimulatorBase* base) { + m_base = base; + Simulation::getCurrent()->setDynamicBoundarySimulator(this); +} + +DynamicBoundarySimulator::~DynamicBoundarySimulator() { + delete StrongCouplingBoundarySolver::getCurrent(); +} + +void DynamicBoundarySimulator::initBoundaryData() { + const std::string& sceneFile = SceneConfiguration::getCurrent()->getSceneFile(); + const Utilities::SceneLoader::Scene& scene = SceneConfiguration::getCurrent()->getScene(); + std::string scene_path = FileSystem::getFilePath(sceneFile); + std::string scene_file_name = FileSystem::getFileName(sceneFile); + // no cache for 2D scenes + // 2D sampling is fast, but storing it would require storing the transformation as well + const bool useCache = m_base->getUseParticleCaching() && !scene.sim2D; + Simulation* sim = Simulation::getCurrent(); + + string cachePath = scene_path + "/Cache"; + + for (unsigned int i = 0; i < scene.boundaryModels.size(); i++) { + string meshFileName = scene.boundaryModels[i]->meshFile; + if (FileSystem::isRelativePath(meshFileName)) + meshFileName = FileSystem::normalizePath(scene_path + "/" + scene.boundaryModels[i]->meshFile); + + // check if mesh file has changed + std::string md5FileName = FileSystem::normalizePath(cachePath + "/" + FileSystem::getFileNameWithExt(meshFileName) + ".md5"); + bool md5 = false; + if (useCache) { + string md5Str = FileSystem::getFileMD5(meshFileName); + if (FileSystem::fileExists(md5FileName)) + md5 = FileSystem::checkMD5(md5Str, md5FileName); + } + + DynamicRigidBody* rb = new DynamicRigidBody(); + rb->setIsAnimated(scene.boundaryModels[i]->isAnimated); + TriangleMesh& geo = rb->getGeometry(); + MeshImport::importMesh(meshFileName, geo, Vector3r::Zero(), Matrix3r::Identity(), scene.boundaryModels[i]->scale); + + std::vector boundaryParticles; + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + // if a samples file is given, use this one + if (scene.boundaryModels[i]->samplesFile != "") { + string particleFileName = scene_path + "/" + scene.boundaryModels[i]->samplesFile; + PartioReaderWriter::readParticles(particleFileName, Vector3r::Zero(), Matrix3r::Identity(), scene.boundaryModels[i]->scale[0], boundaryParticles); + } else // if no samples file is given, sample the surface model + { + // Cache sampling + std::string mesh_base_path = FileSystem::getFilePath(scene.boundaryModels[i]->meshFile); + std::string mesh_file_name = FileSystem::getFileName(scene.boundaryModels[i]->meshFile); + + const string resStr = StringTools::real2String(scene.boundaryModels[i]->scale[0]) + "_" + StringTools::real2String(scene.boundaryModels[i]->scale[1]) + "_" + StringTools::real2String(scene.boundaryModels[i]->scale[2]); + const string modeStr = "_m" + std::to_string(scene.boundaryModels[i]->samplingMode); + const string particleFileName = FileSystem::normalizePath(cachePath + "/" + mesh_file_name + "_sb_" + StringTools::real2String(scene.particleRadius) + "_" + resStr + modeStr + ".bgeo"); + + // check MD5 if cache file is available + bool foundCacheFile = false; + + if (useCache) + foundCacheFile = FileSystem::fileExists(particleFileName); + + if (useCache && foundCacheFile && md5) { + PartioReaderWriter::readParticles(particleFileName, Vector3r::Zero(), Matrix3r::Identity(), 1.0, boundaryParticles); + LOG_INFO << "Loaded cached boundary sampling: " << particleFileName; + } + + if (!useCache || !foundCacheFile || !md5) { + if (!scene.sim2D) { + const auto samplePoissonDisk = [&]() { + LOG_INFO << "Poisson disk surface sampling of " << meshFileName; + START_TIMING("Poisson disk sampling"); + PoissonDiskSampling sampling; + sampling.sampleMesh(geo.numVertices(), geo.getVertices().data(), geo.numFaces(), geo.getFaces().data(), scene.particleRadius, 10, 1, boundaryParticles); + STOP_TIMING_AVG; + }; + const auto sampleRegularTriangle = [&]() { + LOG_INFO << "Regular triangle surface sampling of " << meshFileName; + START_TIMING("Regular triangle sampling"); + RegularTriangleSampling sampling; + sampling.sampleMesh(geo.numVertices(), geo.getVertices().data(), geo.numFaces(), geo.getFaces().data(), 1.5f * scene.particleRadius, boundaryParticles); + STOP_TIMING_AVG; + }; + if (SurfaceSamplingMode::PoissonDisk == scene.boundaryModels[i]->samplingMode) + samplePoissonDisk(); + else if (SurfaceSamplingMode::RegularTriangle == scene.boundaryModels[i]->samplingMode) + sampleRegularTriangle(); + else { + LOG_WARN << "Unknown surface sampling method: " << scene.boundaryModels[i]->samplingMode; + LOG_WARN << "Falling back to:"; + sampleRegularTriangle(); + } + } else { + LOG_INFO << "2D regular sampling of " << meshFileName; + START_TIMING("2D regular sampling"); + RegularSampling2D sampling; + sampling.sampleMesh(Matrix3r::Identity(), Vector3r::Zero(), + geo.numVertices(), geo.getVertices().data(), geo.numFaces(), + geo.getFaces().data(), 1.75f * scene.particleRadius, boundaryParticles); + STOP_TIMING_AVG; + } + + // Cache sampling + if (useCache && (FileSystem::makeDir(cachePath) == 0)) { + LOG_INFO << "Save particle sampling: " << particleFileName; + PartioReaderWriter::writeParticles(particleFileName, (unsigned int)boundaryParticles.size(), boundaryParticles.data(), nullptr, scene.particleRadius); + } + } + } + } + + Matrix3r rot = AngleAxisr(scene.boundaryModels[i]->angle, scene.boundaryModels[i]->axis).toRotationMatrix(); + Quaternionr q(rot); + + rb->initBody(scene.boundaryModels[i]->density, scene.boundaryModels[i]->dynamic, scene.boundaryModels[i]->translation, q, + scene.boundaryModels[i]->scale, scene.boundaryModels[i]->velocity, scene.boundaryModels[i]->friction); + + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + BoundaryModel_Akinci2012* bm = new BoundaryModel_Akinci2012(); + bm->initModel(rb, static_cast(boundaryParticles.size()), &boundaryParticles[0]); + sim->addBoundaryModel(bm); + } else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) { + BoundaryModel_Koschier2017* bm = new BoundaryModel_Koschier2017(); + bm->initModel(rb); + sim->addBoundaryModel(bm); + SPH::TriangleMesh& mesh = rb->getGeometry(); + m_base->initDensityMap(mesh.getVertices(), mesh.getFaces(), scene.boundaryModels[i], md5, false, bm); + } else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) { + BoundaryModel_Bender2019* bm = new BoundaryModel_Bender2019(); + bm->initModel(rb); + sim->addBoundaryModel(bm); + SPH::TriangleMesh& mesh = rb->getGeometry(); + m_base->initVolumeMap(mesh.getVertices(), mesh.getFaces(), scene.boundaryModels[i], md5, false, bm); + } + if (useCache && !md5) + FileSystem::writeMD5File(meshFileName, md5FileName); + rb->updateMeshTransformation(); + + } +} + +void DynamicBoundarySimulator::deferredInit() { + Simulation* sim = Simulation::getCurrent(); + sim->performNeighborhoodSearchSort(); + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + m_base->updateBoundaryParticles(true); + Simulation::getCurrent()->updateBoundaryVolume(); + } else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + m_base->updateDMVelocity(); + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + m_base->updateVMVelocity(); + +#ifdef GPU_NEIGHBORHOOD_SEARCH + // copy the particle data to the GPU + sim->getNeighborhoodSearch()->update_point_sets(); +#endif +} + +void DynamicBoundarySimulator::timeStep() { + Simulation* sim = Simulation::getCurrent(); + updateBoundaryForces(); + const unsigned int nObjects = sim->numberOfBoundaryModels(); + #pragma omp parallel for + for (int i = 0; i < nObjects; i++) { + BoundaryModel* bm = sim->getBoundaryModel(i); + RigidBodyObject* rbo = bm->getRigidBodyObject(); + if (rbo->isDynamic()) { + DynamicRigidBody* drb = dynamic_cast(rbo); + drb->timeStep(); + // Apply damping + drb->setVelocity(drb->getVelocity() * (static_cast(1.0) - m_dampingCoeff)); + drb->setAngularVelocity(drb->getAngularVelocity() * (static_cast(1.0) - m_dampingCoeff)); + } + } + if(sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + m_base->updateBoundaryParticles(false); + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + m_base->updateDMVelocity(); + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + m_base->updateVMVelocity(); +} + +void DynamicBoundarySimulator::updateVelocities() { + Simulation* sim = Simulation::getCurrent(); + updateBoundaryForces(); + const unsigned int nObjects = sim->numberOfBoundaryModels(); + #pragma omp parallel for + for (int i = 0; i < nObjects; i++) { + BoundaryModel* bm = sim->getBoundaryModel(i); + RigidBodyObject* rbo = bm->getRigidBodyObject(); + if (rbo->isDynamic()) { + DynamicRigidBody* drb = dynamic_cast(rbo); + drb->updateVelocity(); + } + } +} + +void DynamicBoundarySimulator::reset() { + Simulation* sim = Simulation::getCurrent(); + for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i++) { + BoundaryModel* bm = sim->getBoundaryModel(i); + if (bm->getRigidBodyObject()->isDynamic() || bm->getRigidBodyObject()->isAnimated()) { + ((DynamicRigidBody*)bm->getRigidBodyObject())->reset(); + } + } + + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) + m_base->updateBoundaryParticles(true); + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) + m_base->updateDMVelocity(); + else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Bender2019) + m_base->updateVMVelocity(); +} diff --git a/Simulator/DynamicBoundarySimulator.h b/Simulator/DynamicBoundarySimulator.h new file mode 100644 index 00000000..7ea3a063 --- /dev/null +++ b/Simulator/DynamicBoundarySimulator.h @@ -0,0 +1,50 @@ +#ifndef __DynamicBoundarySimulator_h__ +#define __DynamicBoundarySimualtor_h__ +#include "BoundarySimulator.h" +namespace SPH { + class SimulatorBase; + class TriangleMesh; + + class DynamicBoundarySimulator : public BoundarySimulator { + private: + Real m_dampingCoeff = 0.0; + + + protected: + SimulatorBase* m_base; + + public: + DynamicBoundarySimulator(SimulatorBase* base); + virtual ~DynamicBoundarySimulator(); + + virtual void initBoundaryData(); + /** This function is called after the simulation scene is loaded and all + * parameters are initialized. While reading a scene file several parameters + * can change. The deferred init function should initialize all values which + * depend on these parameters. + */ + virtual void deferredInit(); + + virtual void timeStep(); + virtual void reset(); + + /** + * This function is used in DFSPH to compute the intermediate velocities after divergence solver + */ + void updateVelocities(); + + FORCE_INLINE const Real& getDampingCoeff() const { + return m_dampingCoeff; + } + FORCE_INLINE Real& getDampingCoeff() { + return m_dampingCoeff; + } + FORCE_INLINE void setDampingCoeff(const Real& value) { + m_dampingCoeff = value; + } + + }; + +} + +#endif diff --git a/Simulator/Exporter/RigidBodyParticleExporter_VTK.cpp b/Simulator/Exporter/RigidBodyParticleExporter_VTK.cpp new file mode 100644 index 00000000..7ef964e6 --- /dev/null +++ b/Simulator/Exporter/RigidBodyParticleExporter_VTK.cpp @@ -0,0 +1,245 @@ +#include "RigidBodyParticleExporter_VTK.h" +#include +#include +#include "SPlisHSPlasH/Simulation.h" +#include +using namespace SPH; +using namespace Utilities; + +RigidBodyParticleExporter_VTK::RigidBodyParticleExporter_VTK(SimulatorBase* base) : ExporterBase(base) { + m_isFirstFrame = true; +} + +RigidBodyParticleExporter_VTK::~RigidBodyParticleExporter_VTK(void) {} + +void RigidBodyParticleExporter_VTK::init(const std::string& outputPath) { + m_exportPath = FileSystem::normalizePath(outputPath + "/vtk"); +} + +void RigidBodyParticleExporter_VTK::step(const unsigned int frame) { + if (!m_active) { + return; + } + writeRigidBodies(frame); +} + +void RigidBodyParticleExporter_VTK::reset() { + m_isFirstFrame = true; +} + +void RigidBodyParticleExporter_VTK::setActive(const bool active) { + ExporterBase::setActive(active); + if (m_active) { + FileSystem::makeDirs(m_exportPath); + } +} + +void RigidBodyParticleExporter_VTK::writeRigidBodies(const unsigned int frame) { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nBoundaryModels = sim->numberOfBoundaryModels(); + + // check if we have a static model + bool isStatic = true; + for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i++) { + BoundaryModel* bm = sim->getBoundaryModel(i); + if (bm->getRigidBodyObject()->isDynamic() || bm->getRigidBodyObject()->isAnimated()) { + isStatic = false; + break; + } + } + +#ifdef USE_DOUBLE + const char* real_str = " double\n"; +#else + const char* real_str = " float\n"; +#endif + + if (m_isFirstFrame || !isStatic) { + for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i++) { + std::string fileName = "rb_particle_data_"; + fileName = fileName + std::to_string(i) + "_" + std::to_string(frame) + ".vtk"; + std::string exportFileName = FileSystem::normalizePath(m_exportPath + "/" + fileName); + + // Open the file + std::ofstream outfile(exportFileName, std::ios::binary); + if (!outfile) { + LOG_WARN << "Cannot open a file to save VTK mesh."; + return; + } + + // Header + outfile << "# vtk DataFile Version 4.2\n"; + outfile << "SPlisHSPlasH rigid body particle data\n"; + outfile << "BINARY\n"; + outfile << "DATASET UNSTRUCTURED_GRID\n"; + + // add attributes + m_attributes.clear(); + StringTools::tokenize(m_base->getValue(SimulatorBase::RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES), m_attributes, ";"); + + ////////////////////////////////////////////////////////////////////////// + // positions and ids exported anyways + m_attributes.erase( + std::remove_if(m_attributes.begin(), m_attributes.end(), [](const std::string& s) { return (s == "position" || s == "id"); }), + m_attributes.end()); + + BoundaryModel_Akinci2012* model = dynamic_cast(sim->getBoundaryModel(i)); + if (model != nullptr) { + Simulation* sim = Simulation::getCurrent(); + const unsigned int nFluids = sim->numberOfFluidModels(); + + const unsigned int numBoundaryParticles = model->numberOfParticles(); + + std::vector positions; + positions.reserve(numBoundaryParticles); + std::vector cells; + cells.reserve(numBoundaryParticles); + std::vector cellTypes; + cellTypes.reserve(numBoundaryParticles); + std::vector attrData; + attrData.reserve(numBoundaryParticles); + + unsigned int counter = 0; + for (unsigned int i = 0; i < numBoundaryParticles; i++) { + + + ////////////////////////////////////////////////////////////////////////// + // export position attribute as POINTS + { + positions.push_back(model->getPosition(i)); + // swap endianess + for (unsigned int c = 0; c < 3; c++) + swapByteOrder(&positions[positions.size() - 1][c]); + } + + ////////////////////////////////////////////////////////////////////////// + // export particle IDs as CELLS + { + unsigned int nodes_per_cell_swapped = 1; + swapByteOrder(&nodes_per_cell_swapped); + unsigned int idSwapped = counter++; + swapByteOrder(&idSwapped); + cells.push_back({ nodes_per_cell_swapped, idSwapped }); + } + ////////////////////////////////////////////////////////////////////////// + // export cell types + { + // the type of a particle cell is always 1 + int cellTypeSwapped = 1; + swapByteOrder(&cellTypeSwapped); + cellTypes.push_back(cellTypeSwapped); + } + + ////////////////////////////////////////////////////////////////////////// + // write additional attributes as per-particle data + //{ + // unsigned int id = model->getParticleId(i); + // swapByteOrder(&id); + // attrData.push_back(id); + //} + } + + //createParticleFile(fileName, model); + + if (outfile.good()) { + // export to vtk + const unsigned int nPoints = (unsigned int)positions.size(); + outfile << "POINTS " << nPoints << real_str; + outfile.write(reinterpret_cast(positions[0].data()), 3 * nPoints * sizeof(Real)); + outfile << "\n"; + + // particles are cells with one element and the index of the particle + outfile << "CELLS " << nPoints << " " << 2 * nPoints << "\n"; + outfile.write(reinterpret_cast(cells[0].data()), 2 * nPoints * sizeof(unsigned int)); + outfile << "\n"; + + outfile << "CELL_TYPES " << nPoints << "\n"; + outfile.write(reinterpret_cast(cellTypes.data()), nPoints * sizeof(int)); + outfile << "\n"; + + outfile << "POINT_DATA " << nPoints << "\n"; + // write IDs + outfile << "SCALARS id unsigned_int 1\n"; + outfile << "LOOKUP_TABLE id_table\n"; + outfile.write(reinterpret_cast(attrData.data()), nPoints * sizeof(unsigned int)); + outfile << "\n"; + + + + ////////////////////////////////////////////////////////////////////////// + // per point fields (all attributes except for positions) + const auto numFields = m_attributes.size(); + outfile << "FIELD FieldData " << std::to_string(numFields) << "\n"; + + // iterate over attributes. + for (const std::string& a : m_attributes) { + const FieldDescription& field = model->getField(a); + + std::string attrNameVTK; + std::regex_replace(std::back_inserter(attrNameVTK), a.begin(), a.end(), std::regex("\\s+"), "_"); + + if (field.type == FieldType::Scalar) { + // write header information + outfile << attrNameVTK << " 1 " << nPoints << real_str; + + // copy data + std::vector attrData; + attrData.reserve(nPoints); + for (unsigned int i = 0u; i < numBoundaryParticles; i++) { + Real val = *((Real*)field.getFct(i)); + swapByteOrder(&val); // swap endianess + attrData.emplace_back(val); + } + + // export to vtk + outfile.write(reinterpret_cast(attrData.data()), nPoints * sizeof(Real)); + } else if (field.type == FieldType::Vector3) { + // write header information + outfile << attrNameVTK << " 3 " << nPoints << real_str; + + // copy from partio data + std::vector attrData; + attrData.reserve(nPoints); + for (unsigned int i = 0u; i < numBoundaryParticles; i++) { + Vector3r val((Real*)field.getFct(i)); + for (unsigned int c = 0; c < 3; c++) + swapByteOrder(&val[c]); // swap endianess + attrData.emplace_back(val); + } + // export to vtk + outfile.write(reinterpret_cast(attrData[0].data()), 3 * nPoints * sizeof(Real)); + } else if (field.type == FieldType::UInt) { + // write header information + outfile << attrNameVTK << " 1 " << nPoints << " unsigned_int\n"; + + // copy data + std::vector attrData; + attrData.reserve(nPoints); + for (unsigned int i = 0u; i < numBoundaryParticles; i++) { + unsigned int val = *((unsigned int*)field.getFct(i)); + swapByteOrder(&val); // swap endianess + attrData.emplace_back(val); + } + // export to vtk + outfile.write(reinterpret_cast(attrData.data()), nPoints * sizeof(unsigned int)); + } + // TODO support other field types + else { + LOG_WARN << "Skipping attribute " << a << ", because it is of unsupported type\n"; + continue; + } + // end of block + outfile << "\n"; + + } + outfile.close(); + } + + } + + } + + m_isFirstFrame = false; + + } +} \ No newline at end of file diff --git a/Simulator/Exporter/RigidBodyParticleExporter_VTK.h b/Simulator/Exporter/RigidBodyParticleExporter_VTK.h new file mode 100644 index 00000000..a523f690 --- /dev/null +++ b/Simulator/Exporter/RigidBodyParticleExporter_VTK.h @@ -0,0 +1,40 @@ +#ifndef __RigidBodyParticleExporter_VTK_h__ +#define __RigidBodyParticleExporter_VTK_h__ + +#include "ExporterBase.h" +#include "SPlisHSPlasH/FluidModel.h" + +namespace SPH { + /** \brief Rigid body particle exporter for the VTK format. + */ + class RigidBodyParticleExporter_VTK : public ExporterBase { + protected: + bool m_isFirstFrame; + std::vector m_attributes; + std::string m_exportPath; + + void writeRigidBodies(const unsigned int frame); + + // VTK expects big endian + template + inline void swapByteOrder(T* v) { + constexpr size_t n = sizeof(T); + uint8_t* bytes = reinterpret_cast(v); + for (unsigned int c = 0u; c < n / 2; c++) + std::swap(bytes[c], bytes[n - c - 1]); + } + + public: + RigidBodyParticleExporter_VTK(SimulatorBase* base); + RigidBodyParticleExporter_VTK(const RigidBodyParticleExporter_VTK&) = delete; + RigidBodyParticleExporter_VTK& operator=(const RigidBodyParticleExporter_VTK&) = delete; + virtual ~RigidBodyParticleExporter_VTK(void); + + virtual void init(const std::string& outputPath); + virtual void step(const unsigned int frame); + virtual void reset(); + virtual void setActive(const bool active); + }; +} + +#endif diff --git a/Simulator/ExporterRegistration.cpp b/Simulator/ExporterRegistration.cpp index 34327646..d726a8e0 100644 --- a/Simulator/ExporterRegistration.cpp +++ b/Simulator/ExporterRegistration.cpp @@ -6,6 +6,7 @@ #include "Exporter/RigidBodyExporter_OBJ.h" #include "Exporter/RigidBodyExporter_PLY.h" #include "Exporter/RigidBodyExporter_VTK.h" +#include "Exporter/RigidBodyParticleExporter_VTK.h" using namespace SPH; @@ -18,4 +19,5 @@ void SimulatorBase::createExporters() addRigidBodyExporter("enableRigidBodyOBJExport", "Rigid Body OBJ Exporter", "Enable/disable rigid body OBJ export.", new RigidBodyExporter_OBJ(this)); addRigidBodyExporter("enableRigidBodyPLYExport", "Rigid Body PLY Exporter", "Enable/disable rigid body PLY export.", new RigidBodyExporter_PLY(this)); addRigidBodyExporter("enableRigidBodyVTKExport", "Rigid Body VTK Exporter", "Enable/disable rigid body VTK export.", new RigidBodyExporter_VTK(this)); + addRigidBodyExporter("enableRigidBodyParticleVTKExport", "Rigid Body Particle VTK Exporter", "Enable/disable rigid body particle VTK export.", new RigidBodyParticleExporter_VTK(this)); } \ No newline at end of file diff --git a/Simulator/GUI/imgui/Simulator_GUI_imgui.cpp b/Simulator/GUI/imgui/Simulator_GUI_imgui.cpp index 293751b0..05dfc050 100644 --- a/Simulator/GUI/imgui/Simulator_GUI_imgui.cpp +++ b/Simulator/GUI/imgui/Simulator_GUI_imgui.cpp @@ -9,16 +9,19 @@ #include "Utilities/FileSystem.h" #include "Simulator/SceneConfiguration.h" #include "LogWindow.h" - +#include "SPlisHSPlasH/DynamicRigidBody.h" #include "imgui.h" #include "imgui_internal.h" #include "backends/imgui_impl_glfw.h" #include "backends/imgui_impl_opengl3.h" +#include "Simulator/DynamicBoundarySimulator.h" +#include "SPlisHSPlasH/StrongCouplingBoundarySolver.h" #define GLFW_INCLUDE_NONE #include + using namespace SPH; using namespace Utilities; @@ -30,6 +33,7 @@ Simulator_GUI_imgui::Simulator_GUI_imgui(SimulatorBase *base) : m_vsync = false; m_iniFound = false; m_showLogWindow = true; + m_currentBoundaryModel = 0; } Simulator_GUI_imgui::~Simulator_GUI_imgui(void) @@ -498,7 +502,113 @@ void Simulator_GUI_imgui::initSimulationParameterGUI() imguiParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getVorticityBase()); imguiParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getElasticityBase()); } - + + if (!m_simulatorBase->isStaticScene() && sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + DynamicBoundarySimulator* simulator = sim->getDynamicBoundarySimulator(); + // Fields for all boundary models + { + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { + imguiParameters::imguiNumericParameter* damping = new imguiParameters::imguiNumericParameter(); + damping->description = "Damping of the system"; + damping->label = "Damping"; + damping->getFct = [simulator]()->Real {return simulator->getDampingCoeff(); }; + damping->setFct = [simulator](Real v) {simulator->setDampingCoeff(v); }; + imguiParameters::addParam("Boundary Model", "General", damping); + + StrongCouplingBoundarySolver* bs = StrongCouplingBoundarySolver::getCurrent(); + + imguiParameters::createParameterObjectGUI(bs); + + imguiParameters::imguiNumericParameter* maxIteration = new imguiParameters::imguiNumericParameter; + maxIteration->description = "max iteration to solve the strong coupling method"; + maxIteration->label = "Max Iterations"; + maxIteration->getFct = [bs]()->unsigned int {return bs->getMaxIterations(); }; + maxIteration->setFct = [bs](unsigned int v) {bs->setMaxIterations(v); }; + imguiParameters::addParam("Boundary Model", "General", maxIteration); + + imguiParameters::imguiNumericParameter* minIteration = new imguiParameters::imguiNumericParameter; + minIteration->description = "min iteration to solve the strong coupling method"; + minIteration->label = "min Iterations"; + minIteration->getFct = [bs]()->unsigned int {return bs->getMinIterations(); }; + minIteration->setFct = [bs](unsigned int v) {bs->setMinIterations(v); }; + imguiParameters::addParam("Boundary Model", "General", minIteration); + + imguiParameters::imguiNumericParameter* maxDensityDeviation = new imguiParameters::imguiNumericParameter; + maxDensityDeviation->description = "max density deviation in % to solve the strong coupling method"; + maxDensityDeviation->label = "max density deviation (%)"; + maxDensityDeviation->getFct = [bs]()->Real {return bs->getMaxDensityDeviation() * 100; }; + maxDensityDeviation->setFct = [bs](Real v) {bs->setMaxDensityDeviation(v / 100); }; + imguiParameters::addParam("Boundary Model", "General", maxDensityDeviation); + } + } + + // Enum for all boundary models + if (sim->numberOfBoundaryModels() > 0) { + BoundaryModel* model = sim->getBoundaryModel(m_currentBoundaryModel); + // Select boundary model + { + imguiParameters::imguiEnumParameter* param = new imguiParameters::imguiEnumParameter(); + param->description = "Select a boundary model to set its parameters below."; + param->label = "Current boundary model"; + param->readOnly = false; + for (unsigned int j = 0; j < sim->numberOfBoundaryModels(); j++) { + param->items.push_back(std::to_string(j)); + } + param->getFct = [this]() -> int { return m_currentBoundaryModel; }; + param->setFct = [this](int v) { m_currentBoundaryModel = v; initSimulationParameterGUI(); }; + imguiParameters::addParam("Boundary Model", "Property", param); + } + + // Show basic properties + { + imguiParameters::imguiBoolParameter* isDynamic = new imguiParameters::imguiBoolParameter(); + isDynamic->description = "Density of the rigid body"; + isDynamic->label = "Dynamic"; + isDynamic->readOnly = true; + isDynamic->getFct = [model]()->Real {return model->getRigidBodyObject()->isDynamic(); }; + imguiParameters::addParam("Boundary Model", "Property", isDynamic); + + if (model->getRigidBodyObject()->isDynamic()) { + DynamicRigidBody* drb = static_cast(model->getRigidBodyObject()); + imguiParameters::imguiNumericParameter* frictionCoeff = new imguiParameters::imguiNumericParameter(); + frictionCoeff->description = "friction coefficient of the rigid body"; + frictionCoeff->label = "friction coefficient"; + frictionCoeff->getFct = [drb]() -> Real {return drb->getFrictionCoeff(); }; + frictionCoeff->setFct = [drb](Real v) {drb->setFrictionCoeff(v); }; + imguiParameters::addParam("Boundary Model", "Property", frictionCoeff); + + imguiParameters::imguiNumericParameter* density = new imguiParameters::imguiNumericParameter(); + density->description = "Density of the rigid body"; + density->label = "Density"; + density->getFct = [drb]() -> Real {return drb->getDensity(); }; + density->setFct = [drb](Real v) {drb->setDensity(v); }; + imguiParameters::addParam("Boundary Model", "Property", density); + + imguiParameters::imguiVec3rParameter* scale = new imguiParameters::imguiVec3rParameter(); + scale->description = "Scale of the rigid body"; + scale->label = "Scale"; + scale->getFct = [drb]()-> Vector3r {return drb->getScale(); }; + scale->readOnly = true; + imguiParameters::addParam("Boundary Model", "Property", scale); + + imguiParameters::imguiVec3rParameter* velocity = new imguiParameters::imguiVec3rParameter(); + velocity->description = "Velocity of the rigid body"; + velocity->label = "Velocity"; + velocity->getFct = [drb]()-> Vector3r {return drb->getVelocity(); }; + velocity->readOnly = true; + imguiParameters::addParam("Boundary Model", "Property", velocity); + + imguiParameters::imguiVec3rParameter* angularVelocity = new imguiParameters::imguiVec3rParameter(); + angularVelocity->description = "Angular Velocity of the rigid body"; + angularVelocity->label = "Angular velocity"; + angularVelocity->getFct = [drb]()-> Vector3r {return drb->getAngularVelocity(); }; + angularVelocity->readOnly = true; + imguiParameters::addParam("Boundary Model", "Property", angularVelocity); + } + } + } + } + initImguiParameters(); if (m_simulatorBase) @@ -639,7 +749,7 @@ void Simulator_GUI_imgui::renderBoundary() const int renderWalls = base->getValue(SimulatorBase::RENDER_WALLS); if (((renderWalls == 1) || (renderWalls == 2)) && - (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012)) + (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019)) { for (int body = sim->numberOfBoundaryModels() - 1; body >= 0; body--) { diff --git a/Simulator/GUI/imgui/Simulator_GUI_imgui.h b/Simulator/GUI/imgui/Simulator_GUI_imgui.h index c5dc7416..be34bed1 100644 --- a/Simulator/GUI/imgui/Simulator_GUI_imgui.h +++ b/Simulator/GUI/imgui/Simulator_GUI_imgui.h @@ -50,6 +50,8 @@ namespace SPH bool m_iniFound; LogWindow* m_logWindow; const float m_baseSize = 15.0f; + // gui for boundary + unsigned int m_currentBoundaryModel; std::vector>& getSelectedParticles() { return m_selectedParticles; } void initImgui(); diff --git a/Simulator/SimulatorBase.cpp b/Simulator/SimulatorBase.cpp index 04109681..afb40b3f 100644 --- a/Simulator/SimulatorBase.cpp +++ b/Simulator/SimulatorBase.cpp @@ -26,6 +26,7 @@ #include "StaticBoundarySimulator.h" #include "Exporter/ExporterBase.h" #include "SPlisHSPlasH/Utilities/MeshImport.h" +#include "DynamicBoundarySimulator.h" #if USE_NFD_FILE_DIALOG #include "extern/nfd/include/nfd.h" #endif @@ -49,6 +50,7 @@ int SimulatorBase::STOP_AT = -1; int SimulatorBase::NUM_STEPS_PER_RENDER = -1; int SimulatorBase::DATA_EXPORT_FPS = -1; int SimulatorBase::PARTICLE_EXPORT_ATTRIBUTES = -1; +int SimulatorBase::RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES = -1; int SimulatorBase::STATE_EXPORT = -1; int SimulatorBase::STATE_EXPORT_FPS = -1; int SimulatorBase::ASYNC_EXPORT = -1; @@ -96,7 +98,8 @@ SimulatorBase::SimulatorBase() m_colorMapType.resize(1, 0); m_renderMinValue.resize(1, 0.0); m_renderMaxValue.resize(1, 5.0); - m_particleAttributes = "velocity"; + m_particleAttributes = "velocity;density;"; + m_rigidbodyParticleAttributes = "body position"; m_timeStepCB = nullptr; m_resetCB = nullptr; m_updateGUI = false; @@ -182,6 +185,11 @@ void SimulatorBase::initParameters() setGroup(PARTICLE_EXPORT_ATTRIBUTES, "Export|General"); setDescription(PARTICLE_EXPORT_ATTRIBUTES, "Attributes that are exported in the partio files (except id and position)."); + RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES = createStringParameter("rigidbodyParticleAttributes", "Rigidbody export attributes", &m_rigidbodyParticleAttributes); + getParameter(RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES)->setReadOnly(true); + setGroup(RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES, "Export|General"); + setDescription(RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES, "Attributes that are exported in the partio files (except id and position)."); + ASYNC_EXPORT = createBoolParameter("enableAsyncExport", "Asynchronous export", &m_enableAsyncExport); setGroup(ASYNC_EXPORT, "Export|General"); setDescription(ASYNC_EXPORT, "Enable/disable asynchronous export of data. The total performance is faster but disable it when measuring the timings of simulation components. \n\n Currently only supported by partio exporter."); @@ -413,6 +421,7 @@ void SimulatorBase::init(int argc, char **argv, const std::string &windowName) ////////////////////////////////////////////////////////////////////////// // init boundary simulation ////////////////////////////////////////////////////////////////////////// + m_boundaryHandlingMethod = scene.boundaryHandlingMethod; m_isStaticScene = true; for (unsigned int i = 0; i < scene.boundaryModels.size(); i++) { @@ -430,7 +439,12 @@ void SimulatorBase::init(int argc, char **argv, const std::string &windowName) else { LOG_INFO << "Initialize dynamic boundary simulation"; - m_boundarySimulator = new PBDBoundarySimulator(this); + if (m_boundaryHandlingMethod == 3) { + // Gissler 2019 + m_boundarySimulator = new DynamicBoundarySimulator(this); + } else { + m_boundarySimulator = new PBDBoundarySimulator(this); + } } initExporters(); @@ -530,6 +544,7 @@ void SimulatorBase::initSimulation() #endif }); } + updateScalarField(); } @@ -538,7 +553,7 @@ void SimulatorBase::deferredInit() m_boundarySimulator->initBoundaryData(); Simulation* sim = Simulation::getCurrent(); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { unsigned int nBoundaryParticles = 0; for (unsigned int i = 0; i < sim->numberOfBoundaryModels(); i++) @@ -584,7 +599,6 @@ void SimulatorBase::deferredInit() void SimulatorBase::runSimulation() { deferredInit(); - if (getStateFile() != "") { // avoid that command line parameter is overwritten @@ -596,7 +610,6 @@ void SimulatorBase::runSimulation() m_doPause = false; } setCommandLineParameter(); - if (!m_useGUI) { const Real stopAt = getValue(SimulatorBase::STOP_AT); @@ -612,8 +625,9 @@ void SimulatorBase::runSimulation() break; } } - else + else { m_gui->run(); + } } void SimulatorBase::cleanup() @@ -834,7 +848,6 @@ void SimulatorBase::reset() #ifdef USE_DEBUG_TOOLS Simulation::getCurrent()->getDebugTools()->reset(); #endif - m_boundarySimulator->reset(); if (m_gui) m_gui->reset(); @@ -905,11 +918,12 @@ void SimulatorBase::timeStep() for (unsigned int i = 0; i < numSteps; i++) { START_TIMING("SimStep"); + Simulation::getCurrent()->getTimeStep()->step(); STOP_TIMING_AVG; - + START_TIMING("SimStepBoundary"); m_boundarySimulator->timeStep(); - + STOP_TIMING_AVG; step(); INCREASE_COUNTER("Time step size", TimeManager::getCurrent()->getTimeStepSize()); @@ -1844,7 +1858,7 @@ void SPH::SimulatorBase::loadState(const std::string &stateFile) } if (dynamic) { - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) updateBoundaryParticles(true); else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) updateDMVelocity(); @@ -2176,7 +2190,7 @@ void SimulatorBase::readFluidParticlesState(const std::string &fileName, FluidMo void SimulatorBase::writeBoundaryState(const std::string &fileName, BoundaryModel *bm) { Simulation *sim = Simulation::getCurrent(); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { BoundaryModel_Akinci2012 *model = static_cast(bm); Partio::ParticlesDataMutable& particleData = *Partio::create(); @@ -2222,7 +2236,7 @@ void SimulatorBase::writeBoundaryState(const std::string &fileName, BoundaryMode void SimulatorBase::readBoundaryState(const std::string &fileName, BoundaryModel *bm) { Simulation *sim = Simulation::getCurrent(); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { if (!FileSystem::fileExists(fileName)) { diff --git a/Simulator/SimulatorBase.h b/Simulator/SimulatorBase.h index c5803f5b..b7802ddc 100644 --- a/Simulator/SimulatorBase.h +++ b/Simulator/SimulatorBase.h @@ -46,6 +46,7 @@ namespace SPH bool m_useParticleCaching; bool m_useGUI; bool m_isStaticScene; + int m_boundaryHandlingMethod; int m_renderWalls; bool m_doPause; Real m_pauseAt; @@ -62,6 +63,7 @@ namespace SPH Vector3r m_cameraPosition; Vector3r m_cameraLookAt; std::string m_particleAttributes; + std::string m_rigidbodyParticleAttributes; std::unique_ptr m_sceneLoader; Real m_nextFrameTime; Real m_nextFrameTimeState; @@ -117,6 +119,7 @@ namespace SPH static int NUM_STEPS_PER_RENDER; static int DATA_EXPORT_FPS; static int PARTICLE_EXPORT_ATTRIBUTES; + static int RIGIDBODY_PARTICLE_EXPORT_ATTRIBUTES; static int STATE_EXPORT; static int STATE_EXPORT_FPS; static int ASYNC_EXPORT; @@ -199,6 +202,7 @@ namespace SPH void setUseParticleCaching(bool val) { m_useParticleCaching = val; } bool getUseGUI() const { return m_useGUI; } void setUseGUI(bool val) { m_useGUI = val; } + int getBoundaryHandlingMethod() const { return m_boundaryHandlingMethod; } const std::string& getColorField(const unsigned int fluidModelIndex) { return m_colorField[fluidModelIndex]; } void setColorField(const unsigned int fluidModelIndex, const std::string& fieldName) { m_colorField[fluidModelIndex] = fieldName; } diff --git a/Simulator/StaticBoundarySimulator.cpp b/Simulator/StaticBoundarySimulator.cpp index 6a1b9cd4..dc69a729 100644 --- a/Simulator/StaticBoundarySimulator.cpp +++ b/Simulator/StaticBoundarySimulator.cpp @@ -61,7 +61,7 @@ void StaticBoundarySimulator::initBoundaryData() MeshImport::importMesh(meshFileName, geo, Vector3r::Zero(), Matrix3r::Identity(), scene.boundaryModels[i]->scale); std::vector boundaryParticles; - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { // if a samples file is given, use this one if (scene.boundaryModels[i]->samplesFile != "") @@ -150,7 +150,7 @@ void StaticBoundarySimulator::initBoundaryData() rb->setRotation0(q); rb->setRotation(q); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { BoundaryModel_Akinci2012 *bm = new BoundaryModel_Akinci2012(); bm->initModel(rb, static_cast(boundaryParticles.size()), &boundaryParticles[0]); @@ -182,7 +182,7 @@ void StaticBoundarySimulator::deferredInit() { Simulation* sim = Simulation::getCurrent(); sim->performNeighborhoodSearchSort(); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) { m_base->updateBoundaryParticles(true); Simulation::getCurrent()->updateBoundaryVolume(); @@ -201,7 +201,7 @@ void StaticBoundarySimulator::deferredInit() void StaticBoundarySimulator::timeStep() { Simulation* sim = Simulation::getCurrent(); - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) m_base->updateBoundaryParticles(false); else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) m_base->updateDMVelocity(); @@ -221,7 +221,7 @@ void StaticBoundarySimulator::reset() } } - if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012) + if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Akinci2012 || sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Gissler2019) m_base->updateBoundaryParticles(true); else if (sim->getBoundaryHandlingMethod() == BoundaryHandlingMethods::Koschier2017) m_base->updateDMVelocity(); diff --git a/Simulator/main.cpp b/Simulator/main.cpp index 8dd5f33e..206b26ec 100644 --- a/Simulator/main.cpp +++ b/Simulator/main.cpp @@ -24,13 +24,15 @@ int main( int argc, char **argv ) REPORT_MEMORY_LEAKS; base = new SimulatorBase(); base->init(argc, argv, "SPlisHSPlasH"); - if (base->getUseGUI()) { if (base->isStaticScene()) gui = new Simulator_GUI_imgui(base); + else if(base->getBoundaryHandlingMethod() == 3) // Gissler 2019 + gui = new Simulator_GUI_imgui(base); else gui = new PBD_Simulator_GUI_imgui(base, ((PBDBoundarySimulator*)base->getBoundarySimulator())->getPBDWrapper()); + base->setGui(gui); } base->run(); diff --git a/Utilities/CMakeLists.txt b/Utilities/CMakeLists.txt index dd9d2199..bf38579e 100644 --- a/Utilities/CMakeLists.txt +++ b/Utilities/CMakeLists.txt @@ -26,6 +26,9 @@ add_library(Utilities SystemInfo.h Timing.h Version.h + VolumeIntegration.cpp + VolumeIntegration.h + ) add_dependencies(Utilities partio zlib) diff --git a/Utilities/VolumeIntegration.cpp b/Utilities/VolumeIntegration.cpp new file mode 100644 index 00000000..97ca8fbb --- /dev/null +++ b/Utilities/VolumeIntegration.cpp @@ -0,0 +1,207 @@ + +#include "VolumeIntegration.h" +#include + +using namespace std; +using namespace Eigen; +using namespace Utilities; + +#define SQR(x) ((x)*(x)) +#define CUBE(x) ((x)*(x)*(x)) + +VolumeIntegration::VolumeIntegration(const unsigned int nVertices, const unsigned int nFaces, Vector3r * const vertices, const unsigned int* indices) + : m_nVertices(nVertices), m_nFaces(nFaces), m_indices(indices), m_face_normals(nFaces), m_weights(nFaces) +{ + // compute center of mass + m_x.setZero(); + for (unsigned int i(0); i < m_nVertices; ++i) + m_x += vertices[i]; + m_x /= (Real)m_nVertices; + + m_vertices.resize(nVertices); + for (unsigned int i(0); i < m_nVertices; ++i) + m_vertices[i] = vertices[i] - m_x; + + for (unsigned int i(0); i < m_nFaces; ++i) + { + const Vector3r &a = m_vertices[m_indices[3 * i]]; + const Vector3r &b = m_vertices[m_indices[3 * i + 1]]; + const Vector3r &c = m_vertices[m_indices[3 * i + 2]]; + + const Vector3r d1 = b - a; + const Vector3r d2 = c - a; + + m_face_normals[i] = d1.cross(d2); + if (m_face_normals[i].isZero(1.e-10)) + m_face_normals[i].setZero(); + else + m_face_normals[i].normalize(); + + m_weights[i] = -m_face_normals[i].dot(a); + } +} + +void VolumeIntegration::compute_inertia_tensor(Real density) +{ + volume_integrals(); + m_volume = static_cast(T0); + + m_mass = static_cast(density * T0); + + /* compute center of mass */ + m_r[0] = static_cast(T1[0] / T0); + m_r[1] = static_cast(T1[1] / T0); + m_r[2] = static_cast(T1[2] / T0); + + /* compute inertia tensor */ + m_theta(0, 0) = static_cast(density * (T2[1] + T2[2])); + m_theta(1, 1) = static_cast(density * (T2[2] + T2[0])); + m_theta(2, 2) = static_cast(density * (T2[0] + T2[1])); + m_theta(0, 1) = m_theta(1, 0) = -density * static_cast(TP[0]); + m_theta(1, 2) = m_theta(2, 1) = -density * static_cast(TP[1]); + m_theta(2, 0) = m_theta(0, 2) = -density * static_cast(TP[2]); + + /* translate inertia tensor to center of mass */ + m_theta(0, 0) -= m_mass * (m_r[1]*m_r[1] + m_r[2]*m_r[2]); + m_theta(1, 1) -= m_mass * (m_r[2]*m_r[2] + m_r[0]*m_r[0]); + m_theta(2, 2) -= m_mass * (m_r[0]*m_r[0] + m_r[1]*m_r[1]); + m_theta(0, 1) = m_theta(1, 0) += m_mass * m_r[0] * m_r[1]; + m_theta(1, 2) = m_theta(2, 1) += m_mass * m_r[1] * m_r[2]; + m_theta(2, 0) = m_theta(0, 2) += m_mass * m_r[2] * m_r[0]; + + m_r += m_x; +} + + +void VolumeIntegration::projection_integrals(unsigned int f) +{ + Real a0, a1, da; + Real b0, b1, db; + Real a0_2, a0_3, a0_4, b0_2, b0_3, b0_4; + Real a1_2, a1_3, b1_2, b1_3; + Real C1, Ca, Caa, Caaa, Cb, Cbb, Cbbb; + Real Cab, Kab, Caab, Kaab, Cabb, Kabb; + + P1 = Pa = Pb = Paa = Pab = Pbb = Paaa = Paab = Pabb = Pbbb = 0.0; + + for (int i = 0; i < 3; i++) + { + a0 = m_vertices[m_indices[3 * f + i]][A]; + b0 = m_vertices[m_indices[3 * f + i]][B]; + a1 = m_vertices[m_indices[3 * f + ((i + 1) % 3)]][A]; + b1 = m_vertices[m_indices[3 * f + ((i + 1) % 3)]][B]; + + da = a1 - a0; + db = b1 - b0; + a0_2 = a0 * a0; a0_3 = a0_2 * a0; a0_4 = a0_3 * a0; + b0_2 = b0 * b0; b0_3 = b0_2 * b0; b0_4 = b0_3 * b0; + a1_2 = a1 * a1; a1_3 = a1_2 * a1; + b1_2 = b1 * b1; b1_3 = b1_2 * b1; + + C1 = a1 + a0; + Ca = a1*C1 + a0_2; Caa = a1*Ca + a0_3; Caaa = a1*Caa + a0_4; + Cb = b1*(b1 + b0) + b0_2; Cbb = b1*Cb + b0_3; Cbbb = b1*Cbb + b0_4; + Cab = 3 * a1_2 + 2 * a1*a0 + a0_2; Kab = a1_2 + 2 * a1*a0 + 3 * a0_2; + Caab = a0*Cab + 4 * a1_3; Kaab = a1*Kab + 4 * a0_3; + Cabb = 4 * b1_3 + 3 * b1_2*b0 + 2 * b1*b0_2 + b0_3; + Kabb = b1_3 + 2 * b1_2*b0 + 3 * b1*b0_2 + 4 * b0_3; + + P1 += db*C1; + Pa += db*Ca; + Paa += db*Caa; + Paaa += db*Caaa; + Pb += da*Cb; + Pbb += da*Cbb; + Pbbb += da*Cbbb; + Pab += db*(b1*Cab + b0*Kab); + Paab += db*(b1*Caab + b0*Kaab); + Pabb += da*(a1*Cabb + a0*Kabb); + } + + P1 /= 2.0; + Pa /= 6.0; + Paa /= 12.0; + Paaa /= 20.0; + Pb /= -6.0; + Pbb /= -12.0; + Pbbb /= -20.0; + Pab /= 24.0; + Paab /= 60.0; + Pabb /= -60.0; +} + +void VolumeIntegration::face_integrals(unsigned int f) +{ + Real w; + Vector3r n; + Real k1, k2, k3, k4; + + projection_integrals(f); + + w = m_weights[f]; + n = m_face_normals[f]; + k1 = (n[C] == 0) ? 0 : 1 / n[C]; + k2 = k1 * k1; k3 = k2 * k1; k4 = k3 * k1; + + Fa = k1 * Pa; + Fb = k1 * Pb; + Fc = -k2 * (n[A]*Pa + n[B]*Pb + w*P1); + + Faa = k1 * Paa; + Fbb = k1 * Pbb; + Fcc = k3 * (SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb + + w*(2*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faaa = k1 * Paaa; + Fbbb = k1 * Pbbb; + Fccc = -k4 * (CUBE(n[A])*Paaa + 3*SQR(n[A])*n[B]*Paab + + 3*n[A]*SQR(n[B])*Pabb + CUBE(n[B])*Pbbb + + 3*w*(SQR(n[A])*Paa + 2*n[A]*n[B]*Pab + SQR(n[B])*Pbb) + + w*w*(3*(n[A]*Pa + n[B]*Pb) + w*P1)); + + Faab = k1 * Paab; + Fbbc = -k2 * (n[A]*Pabb + n[B]*Pbbb + w*Pbb); + Fcca = k3 * (SQR(n[A])*Paaa + 2*n[A]*n[B]*Paab + SQR(n[B])*Pabb + + w*(2*(n[A]*Paa + n[B]*Pab) + w*Pa)); +} + +void VolumeIntegration::volume_integrals() +{ + Real nx, ny, nz; + + T0 = T1[0] = T1[1] = T1[2] + = T2[0] = T2[1] = T2[2] + = TP[0] = TP[1] = TP[2] = 0; + + for (unsigned int i(0); i < m_nFaces; ++i) + { + Vector3r const& n = m_face_normals[i]; + nx = std::abs(n[0]); + ny = std::abs(n[1]); + nz = std::abs(n[2]); + if (nx > ny && nx > nz) + C = 0; + else + C = (ny > nz) ? 1 : 2; + A = (C + 1) % 3; + B = (A + 1) % 3; + + face_integrals(i); + + T0 += n[0] * ((A == 0) ? Fa : ((B == 0) ? Fb : Fc)); + + T1[A] += n[A] * Faa; + T1[B] += n[B] * Fbb; + T1[C] += n[C] * Fcc; + T2[A] += n[A] * Faaa; + T2[B] += n[B] * Fbbb; + T2[C] += n[C] * Fccc; + TP[A] += n[A] * Faab; + TP[B] += n[B] * Fbbc; + TP[C] += n[C] * Fcca; + } + + T1[0] /= 2; T1[1] /= 2; T1[2] /= 2; + T2[0] /= 3; T2[1] /= 3; T2[2] /= 3; + TP[0] /= 2; TP[1] /= 2; TP[2] /= 2; +} diff --git a/Utilities/VolumeIntegration.h b/Utilities/VolumeIntegration.h new file mode 100644 index 00000000..db394125 --- /dev/null +++ b/Utilities/VolumeIntegration.h @@ -0,0 +1,72 @@ + +#ifndef __VOLUME_INTEGRATION_H__ +#define __VOLUME_INTEGRATION_H__ + +#include +#include +#include + +#include + +namespace Utilities +{ + class VolumeIntegration + { + + private: + int A; + int B; + int C; + + // projection integrals + Real P1, Pa, Pb, Paa, Pab, Pbb, Paaa, Paab, Pabb, Pbbb; + // face integrals + Real Fa, Fb, Fc, Faa, Fbb, Fcc, Faaa, Fbbb, Fccc, Faab, Fbbc, Fcca; + // volume integrals + Real T0; + Real T1[3]; + Real T2[3]; + Real TP[3]; + + public: + + VolumeIntegration(const unsigned int nVertices, const unsigned int nFaces, Vector3r * const vertices, const unsigned int* indices); + + /** Compute inertia tensor for given geometry and given density. + */ + void compute_inertia_tensor(Real density); + + /** Return mass of body. */ + Real getMass() const { return m_mass; } + /** Return volume of body. */ + Real getVolume() const { return m_volume; } + /** Return inertia tensor of body. */ + Matrix3r const& getInertia() const { return m_theta; } + /** Return center of mass. */ + Vector3r const& getCenterOfMass() const { return m_r; } + + private: + + void volume_integrals(); + void face_integrals(unsigned int i); + + /** Compute various integrations over projection of face. + */ + void projection_integrals(unsigned int i); + + + std::vector m_face_normals; + std::vector m_weights; + unsigned int m_nVertices; + unsigned int m_nFaces; + std::vector m_vertices; + const unsigned int* m_indices; + + Real m_mass, m_volume; + Vector3r m_r; + Vector3r m_x; + Matrix3r m_theta; + }; +} + +#endif diff --git a/data/Scenes/Cubes.json b/data/Scenes/Cubes.json new file mode 100644 index 00000000..a2562660 --- /dev/null +++ b/data/Scenes/Cubes.json @@ -0,0 +1,93 @@ +{ + "Configuration": + { + "particleRadius": 0.025, + "numberOfStepsPerRenderUpdate": 2, + "density0": 1000, + "simulationMethod": 0, + "gravitation": [0,-9.81,0], + "cflMethod": 1, + "cflFactor": 0.5, + "cflMaxTimeStepSize": 0.005, + "maxIterations": 100, + "maxError": 0.05, + "maxIterationsV": 100, + "maxErrorV": 0.1, + "stiffness": 50000, + "exponent": 7, + "velocityUpdateMethod": 0, + "enableDivergenceSolver": true, + "boundaryHandlingMethod": 3 + }, + "Simulation": + { + "timeStepSize": 0.005, + "maxIter" : 5, + "maxIterVel" : 5, + "velocityUpdateMethod" : 0, + "contactTolerance": 0.06, + "tetModelSimulationMethod": 2, + "triangleModelSimulationMethod": 2, + "triangleModelBendingMethod": 2, + "contactStiffnessRigidBody" : 1.0, + "contactStiffnessParticleRigidBody": 100.0, + "cloth_stiffness": 1.0, + "cloth_bendingStiffness": 0.005, + "cloth_xxStiffness": 1.0, + "cloth_yyStiffness": 1.0, + "cloth_xyStiffness": 1.0, + "cloth_xyPoissonRatio": 0.3, + "cloth_yxPoissonRatio": 0.3, + "cloth_normalizeStretch": 0, + "cloth_normalizeShear": 0, + "solid_stiffness": 1.0, + "solid_poissonRatio": 0.2, + "solid_normalizeStretch": 0, + "solid_normalizeShear": 0 + }, + "RigidBodies": [ + { + "id": 1, + "geometryFile": "../models/UnitBox.obj", + "translation": [0,4,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0, + "scale": [3.1, 8, 1.6], + "color": [0.1, 0.4, 0.6, 1.0], + "isDynamic": false, + "isWall": true, + "restitution" : 0.6, + "friction" : 0.0, + "collisionObjectType": 2, + "collisionObjectScale": [3.1, 8, 1.6], + "invertSDF": true, + "mapInvert": true, + "mapThickness": 0.0, + "mapResolution": [30,40,15], + "samplingMode": 1 + }, + { + "id": 2, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 400, + "translation": [0,2,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0, + "scale": [0.4, 0.4, 0.4], + "velocity": [0,0,0], + "restitution" : 0.6, + "friction" : 0.2, + "color": [0.3, 0.5, 0.8, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.4, 0.4, 0.4], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20], + "samplingMode": 1 + } + ] +} + + + diff --git a/data/Scenes/CubesTwo.json b/data/Scenes/CubesTwo.json new file mode 100644 index 00000000..84328199 --- /dev/null +++ b/data/Scenes/CubesTwo.json @@ -0,0 +1,113 @@ +{ + "Configuration": + { + "particleRadius": 0.025, + "numberOfStepsPerRenderUpdate": 2, + "density0": 1000, + "simulationMethod": 0, + "gravitation": [0,-9.81,0], + "cflMethod": 1, + "cflFactor": 0.5, + "cflMaxTimeStepSize": 0.005, + "maxIterations": 100, + "maxError": 0.05, + "maxIterationsV": 100, + "maxErrorV": 0.1, + "stiffness": 50000, + "exponent": 7, + "velocityUpdateMethod": 0, + "enableDivergenceSolver": true, + "boundaryHandlingMethod": 3 + }, + "Simulation": + { + "timeStepSize": 0.005, + "maxIter" : 5, + "maxIterVel" : 5, + "velocityUpdateMethod" : 0, + "contactTolerance": 0.06, + "tetModelSimulationMethod": 2, + "triangleModelSimulationMethod": 2, + "triangleModelBendingMethod": 2, + "contactStiffnessRigidBody" : 1.0, + "contactStiffnessParticleRigidBody": 100.0, + "cloth_stiffness": 1.0, + "cloth_bendingStiffness": 0.005, + "cloth_xxStiffness": 1.0, + "cloth_yyStiffness": 1.0, + "cloth_xyStiffness": 1.0, + "cloth_xyPoissonRatio": 0.3, + "cloth_yxPoissonRatio": 0.3, + "cloth_normalizeStretch": 0, + "cloth_normalizeShear": 0, + "solid_stiffness": 1.0, + "solid_poissonRatio": 0.2, + "solid_normalizeStretch": 0, + "solid_normalizeShear": 0 + }, + "RigidBodies": [ + { + "id": 1, + "geometryFile": "../models/UnitBox.obj", + "translation": [0,4,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0, + "scale": [3.1, 8, 1.6], + "color": [0.1, 0.4, 0.6, 1.0], + "isDynamic": false, + "isWall": true, + "restitution" : 0.6, + "friction" : 0.0, + "collisionObjectType": 2, + "collisionObjectScale": [3.1, 8, 1.6], + "invertSDF": true, + "mapInvert": true, + "mapThickness": 0.0, + "mapResolution": [30,40,15], + "samplingMode": 1 + }, + { + "id": 2, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 400, + "translation": [-1,2.0,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0.0, + "scale": [0.2, 0.4, 0.4], + "velocity": [3,0,0], + "restitution" : 0.6, + "friction" : 0.2, + "color": [0.3, 0.5, 0.8, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.4, 0.4, 0.4], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20], + "samplingMode": 1 + }, + { + "id": 3, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 400, + "translation": [1,2.0,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0.0, + "scale": [0.1, 0.1, 0.1], + "velocity": [-3,0,0], + "restitution" : 0.6, + "friction" : 0.2, + "color": [0.3, 0.5, 0.8, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.4, 0.4, 0.4], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20], + "samplingMode": 1 + } + ] +} + + + diff --git a/data/Scenes/DamBreakModel.json b/data/Scenes/DamBreakModel.json index 24e369b4..0529426d 100644 --- a/data/Scenes/DamBreakModel.json +++ b/data/Scenes/DamBreakModel.json @@ -18,7 +18,7 @@ "exponent": 7, "velocityUpdateMethod": 0, "enableDivergenceSolver": true, - "particleAttributes": "density;velocity", + "particleAttributes": "density;velocity;pressure", "boundaryHandlingMethod": 2 }, "Materials": [ diff --git a/data/Scenes/DamBreakWithCubes.json b/data/Scenes/DamBreakWithCubes.json new file mode 100644 index 00000000..7236f8ca --- /dev/null +++ b/data/Scenes/DamBreakWithCubes.json @@ -0,0 +1,136 @@ +{ + "Configuration": + { + "particleRadius": 0.025, + "numberOfStepsPerRenderUpdate": 2, + "density0": 1000, + "simulationMethod": 3, + "gravitation": [0,-9.81,0], + "cflMethod": 1, + "cflFactor": 0.5, + "cflMaxTimeStepSize": 0.005, + "maxIterations": 100, + "maxError": 0.05, + "maxIterationsV": 100, + "maxErrorV": 0.1, + "stiffness": 50000, + "exponent": 7, + "velocityUpdateMethod": 0, + "enableDivergenceSolver": true, + "boundaryHandlingMethod": 3 + }, + "Simulation": + { + "timeStepSize": 0.005, + "maxIter" : 5, + "maxIterVel" : 5, + "velocityUpdateMethod" : 0, + "contactTolerance": 0.06, + "tetModelSimulationMethod": 2, + "triangleModelSimulationMethod": 2, + "triangleModelBendingMethod": 2, + "contactStiffnessRigidBody" : 1.0, + "contactStiffnessParticleRigidBody": 100.0, + "cloth_stiffness": 1.0, + "cloth_bendingStiffness": 0.005, + "cloth_xxStiffness": 1.0, + "cloth_yyStiffness": 1.0, + "cloth_xyStiffness": 1.0, + "cloth_xyPoissonRatio": 0.3, + "cloth_yxPoissonRatio": 0.3, + "cloth_normalizeStretch": 0, + "cloth_normalizeShear": 0, + "solid_stiffness": 1.0, + "solid_poissonRatio": 0.2, + "solid_normalizeStretch": 0, + "solid_normalizeShear": 0 + }, + "RigidBodies": [ + { + "id": 1, + "geometryFile": "../models/UnitBox.obj", + "translation": [0,4,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0, + "scale": [3.1, 8, 1.6], + "color": [0.1, 0.4, 0.6, 1.0], + "isDynamic": false, + "isWall": true, + "restitution" : 0.6, + "friction" : 0.0, + "collisionObjectType": 2, + "collisionObjectScale": [3.1, 8, 1.6], + "invertSDF": true, + "mapInvert": true, + "mapThickness": 0.0, + "mapResolution": [30,40,15] + }, + { + "id": 2, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 400, + "translation": [0.1,2.0,0], + "rotationAxis": [1, 0, 0], + "rotationAngle": 0.0, + "scale": [0.4, 0.4, 0.4], + "velocity": [0,0,0], + "restitution" : 0.6, + "friction" : 0.2, + "color": [0.3, 0.5, 0.8, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.4, 0.4, 0.4], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20] + }, + { + "id": 3, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 700, + "translation": [-0.2, 3, 0], + "rotationAxis": [1, 1, 1], + "rotationAngle": 0.5, + "scale": [0.1, 0.5, 0.1], + "restitution" : 0.6, + "friction" : 0.2, + "color": [1.0, 1.0, 1.0, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.1, 0.5, 0.1], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20] + }, + { + "id": 4, + "geometryFile": "../models/UnitBox.obj", + "isDynamic": 1, + "density": 100, + "translation": [0.4,1.0,0.0], + "rotationAxis": [0, 0, 1], + "rotationAngle": 1.0, + "scale": [0.2, 0.2, 0.2], + "restitution" : 0.6, + "friction" : 0.2, + "color": [0.5, 0.2, 0.6, 1.0], + "collisionObjectType": 2, + "collisionObjectScale": [0.2, 0.2, 0.2], + "mapInvert": false, + "mapThickness": 0.0, + "mapResolution": [20,20,20] + } + ], + "FluidBlocks": [ + { + "denseMode": 0, + "start": [-0.75, 0.0, -0.75], + "end": [0.75, 1.5, 0.75], + "translation": [-0.75, 0.0, 0.0], + "scale": [1,1,1] + } + ] +} + + + diff --git a/data/Scenes/DamBreakWithObjects.json b/data/Scenes/DamBreakWithObjects.json index b0f507f3..d87712bd 100644 --- a/data/Scenes/DamBreakWithObjects.json +++ b/data/Scenes/DamBreakWithObjects.json @@ -17,7 +17,7 @@ "exponent": 7, "velocityUpdateMethod": 0, "enableDivergenceSolver": true, - "boundaryHandlingMethod": 2 + "boundaryHandlingMethod": 0 }, "Simulation": { diff --git a/data/Scenes/DoubleDamBreakWithSphere.json b/data/Scenes/DoubleDamBreakWithSphere.json index eb4b6ad3..fb6bc07f 100644 --- a/data/Scenes/DoubleDamBreakWithSphere.json +++ b/data/Scenes/DoubleDamBreakWithSphere.json @@ -18,7 +18,7 @@ "exponent": 7, "velocityUpdateMethod": 0, "enableDivergenceSolver": true, - "boundaryHandlingMethod": 2 + "boundaryHandlingMethod": 0 }, "Materials": [ { diff --git a/pySPlisHSPlasH/RigidBodyModule.cpp b/pySPlisHSPlasH/RigidBodyModule.cpp index 9d892218..1d8da191 100644 --- a/pySPlisHSPlasH/RigidBodyModule.cpp +++ b/pySPlisHSPlasH/RigidBodyModule.cpp @@ -5,6 +5,7 @@ #include #include +#include #include #include @@ -64,4 +65,25 @@ void RigidBodyModule(py::module m_sub){ return py::memoryview::from_buffer(base_ptr, { faces.size() }, { sizeof(unsigned int) }); }); + // --------------------------------------- + // Class Dynamic Rigid Body + // --------------------------------------- + py::class_(m_sub, "DynamicRigidBody") + .def(py::init<>()) + .def("setWorldSpacePosition", &SPH::DynamicRigidBody::setWorldSpacePosition) + .def("setWorldSpaceRotation", &SPH::DynamicRigidBody::setWorldSpaceRotation) + .def("animate", &SPH::DynamicRigidBody::animate) + .def("getGeometry", &SPH::DynamicRigidBody::getGeometry) + .def("getVertexBuffer", [](SPH::DynamicRigidBody& obj) -> py::memoryview { + auto vertices = obj.getVertices(); + void* base_ptr = &vertices[0][0]; + int num_vert = vertices.size(); + return py::memoryview::from_buffer((Real*)base_ptr, { num_vert, 3 }, { sizeof(Real) * 3, sizeof(Real) }); + }) + .def("getFaceBuffer", [](SPH::DynamicRigidBody& obj) -> py::memoryview { + auto faces = obj.getFaces(); + unsigned int* base_ptr = faces.data(); + return py::memoryview::from_buffer(base_ptr, { faces.size() }, { sizeof(unsigned int) }); + }); + } diff --git a/pySPlisHSPlasH/UtilitiesModule.cpp b/pySPlisHSPlasH/UtilitiesModule.cpp index 7584ed2a..f6479b84 100644 --- a/pySPlisHSPlasH/UtilitiesModule.cpp +++ b/pySPlisHSPlasH/UtilitiesModule.cpp @@ -277,20 +277,20 @@ void UtilitiesModule(py::module m) { .def(py::init<>()) .def(py::init, std::string, bool, Real, - Eigen::Matrix, unsigned int, bool>(), + Eigen::Matrix, unsigned int, bool, Vector3r, Real, Real>(), "samplesFile"_a = "", "meshFile"_a = "", "translation"_a = Vector3r::Zero(), - "axis"_a = Vector3r(1,0,0), "angle"_a = 0.0, "scale"_a = Vector3r::Ones(), + "axis"_a = Vector3r(1,0,0), "angle"_a = 0.0, "scale"_a = Vector3r::Ones(), "isDynamic"_a = false, "isWall"_a = false, "color"_a = Eigen::Vector4f(1.f, 0.f, 0.f, 0.f), "mapFile"_a = "", "mapInvert"_a = false, "mapThickness"_a = 0.0, "mapResolution"_a = Eigen::Matrix(20, 20, 20), - "samplingMode"_a = 0, "isAnimated"_a = false) + "samplingMode"_a = 0, "isAnimated"_a = false, "velocity"_a = Vector3r::Zero(), "density"_a = 1000, "friction"_a = 0.2) .def_readwrite("samplesFile", &Utilities::BoundaryParameterObject::samplesFile) .def_readwrite("meshFile", &Utilities::BoundaryParameterObject::meshFile) .def_readwrite("translation", &Utilities::BoundaryParameterObject::translation) - .def_readwrite("axis", &Utilities::BoundaryParameterObject::axis) + .def_readwrite("axis", &Utilities::BoundaryParameterObject::axis) .def_readwrite("angle", &Utilities::BoundaryParameterObject::angle) - .def_readwrite("scale", &Utilities::BoundaryParameterObject::scale) + .def_readwrite("scale", &Utilities::BoundaryParameterObject::scale) .def_readwrite("dynamic", &Utilities::BoundaryParameterObject::dynamic) .def_readwrite("isWall", &Utilities::BoundaryParameterObject::isWall) .def_readwrite("color", &Utilities::BoundaryParameterObject::color) @@ -299,7 +299,10 @@ void UtilitiesModule(py::module m) { .def_readwrite("mapThickness", &Utilities::BoundaryParameterObject::mapThickness) .def_readwrite("mapResolution", &Utilities::BoundaryParameterObject::mapResolution) .def_readwrite("samplingMode", &Utilities::BoundaryParameterObject::samplingMode) - .def_readwrite("isAnimated", &Utilities::BoundaryParameterObject::isAnimated); + .def_readwrite("isAnimated", &Utilities::BoundaryParameterObject::isAnimated) + .def_readwrite("velocity", &Utilities::BoundaryParameterObject::velocity) + .def_readwrite("density", &Utilities::BoundaryParameterObject::density) + .def_readwrite("friction", &Utilities::BoundaryParameterObject::friction); py::class_(m_sub_sub, "FluidData") .def(py::init<>())