From 0d1656a0086a99205f43e204a988c52ef197cde5 Mon Sep 17 00:00:00 2001
From: Jan Bender <bender@cs.rwth-aachen.de>
Date: Fri, 8 Feb 2019 09:40:26 +0100
Subject: [PATCH] - added the Implicit SPH Formulation for Incompressible
Linearly Elastic Solids of Peer et al. 2017 - added Corotated SPH for
deformable solids of Becker et al. 2009 - improved particle coloring
---
CMakeLists.txt | 8 +-
Changelog.txt | 6 +
README.md | 18 +-
SPlisHSPlasH/BoundaryModel.cpp | 8 +-
SPlisHSPlasH/BoundaryModel.h | 18 +-
SPlisHSPlasH/CMakeLists.txt | 21 +-
SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp | 35 +-
SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp | 2 +-
SPlisHSPlasH/Drag/DragForce_Macklin2014.cpp | 2 +-
SPlisHSPlasH/Elasticity/ElasticityBase.cpp | 40 ++
SPlisHSPlasH/Elasticity/ElasticityBase.h | 29 +
.../Elasticity/Elasticity_Becker2009.cpp | 335 +++++++++++
.../Elasticity/Elasticity_Becker2009.h | 59 ++
.../Elasticity/Elasticity_Peer2018.cpp | 562 ++++++++++++++++++
SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h | 73 +++
SPlisHSPlasH/FluidModel.cpp | 106 +++-
SPlisHSPlasH/FluidModel.h | 33 +-
SPlisHSPlasH/IISPH/TimeStepIISPH.cpp | 31 +-
SPlisHSPlasH/PBF/SimulationDataPBF.cpp | 4 +-
SPlisHSPlasH/PBF/TimeStepPBF.cpp | 21 +-
SPlisHSPlasH/PCISPH/SimulationDataPCISPH.cpp | 2 +-
SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp | 23 +-
SPlisHSPlasH/Simulation.cpp | 23 +-
SPlisHSPlasH/Simulation.h | 2 +-
.../SurfaceTension_Akinci2013.cpp | 7 +-
.../SurfaceTension_Becker2007.cpp | 5 +-
.../SurfaceTension/SurfaceTension_He2014.cpp | 12 +-
SPlisHSPlasH/Utilities/MathFunctions.cpp | 231 +++++++
SPlisHSPlasH/Utilities/MathFunctions.h | 26 +
SPlisHSPlasH/Utilities/SceneLoader.cpp | 29 +
SPlisHSPlasH/Utilities/SceneLoader.h | 10 +
.../Viscosity/Viscosity_Bender2017.cpp | 11 +-
SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp | 6 +-
SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp | 8 +-
SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp | 7 +-
.../Viscosity/Viscosity_Takahashi2015.cpp | 10 +-
.../Viscosity/Viscosity_Weiler2018.cpp | 31 +-
SPlisHSPlasH/Viscosity/Viscosity_XSPH.cpp | 2 +-
.../Vorticity/MicropolarModel_Bender2017.cpp | 66 +-
.../Vorticity/VorticityConfinement.cpp | 4 +
SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp | 19 +-
{Demos => Simulators}/CMakeLists.txt | 2 +-
.../Common/SimulatorBase.cpp | 397 ++++++++-----
.../Common/SimulatorBase.h | 48 +-
.../Common/TweakBarParameters.cpp | 45 +-
.../Common/TweakBarParameters.h | 0
.../DynamicBoundarySimulator}/CMakeLists.txt | 30 +-
.../PBDRigidBody.h | 0
.../PBDWrapper.cpp | 4 +-
.../PositionBasedDynamicsWrapper/PBDWrapper.h | 0
.../DynamicBoundarySimulator}/main.cpp | 132 +++-
.../StaticBoundarySimulator}/CMakeLists.txt | 26 +-
.../StaticBoundarySimulator}/main.cpp | 130 +++-
Visualization/MiniGL.cpp | 31 +-
Visualization/MiniGL.h | 3 +-
bin/Buckling.bat | 2 -
bin/Buckling_Bender2017.bat | 2 +-
bin/Buckling_Peer2015.bat | 2 +-
bin/Buckling_Peer2016.bat | 2 +-
bin/Buckling_Takahashi2015.bat | 2 +-
bin/Buckling_Weiler2018.bat | 2 +-
bin/Bunny_vs_Dragon.bat | 2 +-
bin/Coiling_Bender2017.bat | 2 +-
bin/Coiling_Peer2015.bat | 2 +-
bin/Coiling_Peer2016.bat | 2 +-
bin/Coiling_Takahashi2015.bat | 2 +-
bin/Coiling_Weiler2018.bat | 2 +-
bin/DamBreakModel.bat | 2 +-
bin/DamBreakModelDragons.bat | 2 +-
bin/DamBreakModel_2D.bat | 2 +-
bin/DamBreakWithObjects.bat | 2 +-
bin/DoubleDamBreak.bat | 2 +-
bin/DoubleDamBreakMultiPhase.bat | 2 +
bin/DoubleDamBreakTwoPhase.bat | 2 -
bin/DoubleDamBreakWithSphere.bat | 2 +-
bin/DragTest.bat | 2 +-
bin/Emitter.bat | 2 +-
bin/MultiPhaseColoring.bat | 2 +
bin/Obstacle.bat | 2 +-
bin/ViscousBunny.bat | 2 +-
data/Scenes/DamBreakModel.json | 3 +-
data/Scenes/DeformableModel.json | 55 ++
...ase.json => DoubleDamBreakMultiPhase.json} | 0
data/Scenes/MultiPhaseColoring.json | 124 ++++
data/Scenes/Obstacle.json | 2 +-
85 files changed, 2641 insertions(+), 386 deletions(-)
create mode 100644 SPlisHSPlasH/Elasticity/ElasticityBase.cpp
create mode 100644 SPlisHSPlasH/Elasticity/ElasticityBase.h
create mode 100644 SPlisHSPlasH/Elasticity/Elasticity_Becker2009.cpp
create mode 100644 SPlisHSPlasH/Elasticity/Elasticity_Becker2009.h
create mode 100644 SPlisHSPlasH/Elasticity/Elasticity_Peer2018.cpp
create mode 100644 SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h
create mode 100644 SPlisHSPlasH/Utilities/MathFunctions.cpp
create mode 100644 SPlisHSPlasH/Utilities/MathFunctions.h
rename {Demos => Simulators}/CMakeLists.txt (63%)
rename Demos/Common/DemoBase.cpp => Simulators/Common/SimulatorBase.cpp (70%)
rename Demos/Common/DemoBase.h => Simulators/Common/SimulatorBase.h (68%)
rename {Demos => Simulators}/Common/TweakBarParameters.cpp (89%)
rename {Demos => Simulators}/Common/TweakBarParameters.h (100%)
rename {Demos/DynamicBoundaryDemo => Simulators/DynamicBoundarySimulator}/CMakeLists.txt (75%)
rename {Demos/DynamicBoundaryDemo => Simulators/DynamicBoundarySimulator}/PositionBasedDynamicsWrapper/PBDRigidBody.h (100%)
rename {Demos/DynamicBoundaryDemo => Simulators/DynamicBoundarySimulator}/PositionBasedDynamicsWrapper/PBDWrapper.cpp (99%)
rename {Demos/DynamicBoundaryDemo => Simulators/DynamicBoundarySimulator}/PositionBasedDynamicsWrapper/PBDWrapper.h (100%)
rename {Demos/DynamicBoundaryDemo => Simulators/DynamicBoundarySimulator}/main.cpp (76%)
rename {Demos/StaticBoundaryDemo => Simulators/StaticBoundarySimulator}/CMakeLists.txt (66%)
rename {Demos/StaticBoundaryDemo => Simulators/StaticBoundarySimulator}/main.cpp (74%)
delete mode 100644 bin/Buckling.bat
create mode 100644 bin/DoubleDamBreakMultiPhase.bat
delete mode 100644 bin/DoubleDamBreakTwoPhase.bat
create mode 100644 bin/MultiPhaseColoring.bat
create mode 100644 data/Scenes/DeformableModel.json
rename data/Scenes/{DoubleDamBreakTwoPhase.json => DoubleDamBreakMultiPhase.json} (100%)
create mode 100644 data/Scenes/MultiPhaseColoring.json
diff --git a/CMakeLists.txt b/CMakeLists.txt
index d37f379..16880c1 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -37,7 +37,7 @@ if (NOT SPH_LIBS_ONLY)
endif()
add_subdirectory(extern/AntTweakBar)
add_subdirectory(extern/glew)
- add_subdirectory(Demos)
+ add_subdirectory(Simulators)
add_subdirectory(Tools)
add_subdirectory(Tests)
endif()
@@ -54,7 +54,7 @@ ExternalProject_Add(
Ext_PBD
PREFIX "${CMAKE_SOURCE_DIR}/extern/PositionBasedDynamics"
GIT_REPOSITORY https://github.com/InteractiveComputerGraphics/PositionBasedDynamics.git
- GIT_TAG "5980708692418f45a26bbd9d3318235a41312e75"
+ GIT_TAG "46f70c2e88ec8c8d7acc15e51c967ed35754f393"
INSTALL_DIR ${ExternalInstallDir}/PositionBasedDynamics
CMAKE_ARGS -DCMAKE_BUILD_TYPE=${EXT_CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX:PATH=${ExternalInstallDir}/PositionBasedDynamics -DPBD_NO_DEMOS:BOOL=1 -DPBD_EXTERNALINSTALLDIR:PATH=${ExternalInstallDir} -DUSE_DOUBLE_PRECISION:BOOL=${USE_DOUBLE_PRECISION}
)
@@ -64,7 +64,7 @@ ExternalProject_Add(
Ext_CompactNSearch
PREFIX "${CMAKE_SOURCE_DIR}/extern/CompactNSearch"
GIT_REPOSITORY https://github.com/InteractiveComputerGraphics/CompactNSearch.git
- GIT_TAG "4ad505f5457e34ec479994740b75767241cff468"
+ GIT_TAG "93724411993689ea4b0fa2daed0b295ebe1e329f"
INSTALL_DIR ${ExternalInstallDir}/CompactNSearch
CMAKE_ARGS -DCMAKE_BUILD_TYPE=${EXT_CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX:PATH=${ExternalInstallDir}/CompactNSearch -DUSE_DOUBLE_PRECISION:BOOL=${USE_DOUBLE_PRECISION}
)
@@ -74,7 +74,7 @@ ExternalProject_Add(
Ext_GenericParameters
PREFIX "${CMAKE_SOURCE_DIR}/extern/GenericParameters"
GIT_REPOSITORY https://github.com/InteractiveComputerGraphics/GenericParameters.git
- GIT_TAG "1ec904bf8555e78ae0d8ba2f9f9a395371c5d4eb"
+ GIT_TAG "b1ad669fac8d106515f6aa8514a03598d5766a36"
INSTALL_DIR ${ExternalInstallDir}/GenericParameters
CMAKE_ARGS -DCMAKE_BUILD_TYPE=${EXT_CMAKE_BUILD_TYPE} -DCMAKE_INSTALL_PREFIX:PATH=${ExternalInstallDir}/GenericParameters -DGENERICPARAMETERS_NO_TESTS:BOOL=1
)
diff --git a/Changelog.txt b/Changelog.txt
index 67e16c8..a9ad7af 100644
--- a/Changelog.txt
+++ b/Changelog.txt
@@ -1,4 +1,10 @@
+2.3.0
+ - added the Implicit SPH Formulation for Incompressible Linearly Elastic Solids of Peer et al. 2017
+ - added Corotated SPH for deformable solids of Becker et al. 2009
- SPlisHSPlasH now supports 2D simulations
+ - SPlisHSPlasH now has enhanced particle coloring
+ - partio export of arbitrary particle attributes is now supported
+ - renamed Static/DynamicBoundaryDemo to Static/DynamicBoundarySimulator
- added documentation of file format
- added colormaps
- fixed race condition
diff --git a/README.md b/README.md
index 712c6d6..f98b401 100644
--- a/README.md
+++ b/README.md
@@ -34,17 +34,21 @@ Note: Please use a 64-bit target on a 64-bit operating system. 32-bit builds on
## Features
SPlisHSPlasH implements:
-* an open-source SPH fluid simulation
+* an open-source SPH fluid simulation (2D & 3D)
* several implicit pressure solvers (WCSPH, PCISPH, PBF, IISPH, DFSPH, PF)
* explicit and implicit viscosity methods
* current surface tension approaches
* different vorticity methods
+* computation of drag forces
* support for multi-phase simulations
+* simulation of deformable solids
+* rigid-fluid coupling with static and dynamic bodies
+* two-way coupling with deformable solids
* fluid emitters
* a json-based scene file importer
* automatic surface sampling
* a tool for volume sampling of closed geometries
-* partio file export
+* partio file export of all particle data
## Pressure Solvers
@@ -102,6 +106,12 @@ The SPlisHSPlasH library implements the drag force computation of the following
* Miles Macklin, Matthias Müller, Nuttapong Chentanez and Tae-Yong Kim. Unified Particle Physics for Real-Time Applications. ACM Trans. Graph., 33(4), 2014
+## Elastic Forces
+
+* M. Becker, M. Ihmsen, and M. Teschner. Corotated SPH for deformable solids. Proceedings of Eurographics Conference on Natural Phenomena, 2009
+
+* A. Peer, C. Gissler, S. Band, and M. Teschner. An Implicit SPH Formulation for Incompressible Linearly Elastic Solids. Computer Graphics Forum, 2017
+
## Multi-Phase Fluid Simulation
The SPlisHSPlasH library implements the following publication to realize multi-phase simulations:
@@ -132,6 +142,8 @@ The following videos were generated using the SPlisHSPlasH library:
* Markus Becker and Matthias Teschner. Weakly compressible SPH for free surface flows. In Proceedings of ACM SIGGRAPH/Eurographics Symposium on Computer Animation, 2007. Eurographics Association.
+* M. Becker, M. Ihmsen, and M. Teschner. Corotated SPH for deformable solids. Proceedings of Eurographics Conference on Natural Phenomena, 2009
+
* Jan Bender and Dan Koschier. Divergence-free smoothed particle hydrodynamics. In Proceedings of ACM SIGGRAPH / Eurographics Symposium on Computer Animation, 2015. ACM.
* Jan Bender and Dan Koschier. Divergence-free SPH for incompressible and viscous fluids. IEEE Transactions on Visualization and Computer Graphics, 2017.
@@ -154,6 +166,8 @@ The following videos were generated using the SPlisHSPlasH library:
* Miles Macklin, Matthias Müller, Nuttapong Chentanez and Tae-Yong Kim. Unified Particle Physics for Real-Time Applications. ACM Trans. Graph., 33(4), 2014
+* A. Peer, C. Gissler, S. Band, and M. Teschner. An Implicit SPH Formulation for Incompressible Linearly Elastic Solids. Computer Graphics Forum, 2017
+
* Andreas Peer, Markus Ihmsen, Jens Cornelis, and Matthias Teschner. An Implicit Viscosity Formulation for SPH Fluids. ACM Trans. Graph., 34(4), 2015.
* Andreas Peer and Matthias Teschner. Prescribed Velocity Gradients for Highly Viscous SPH Fluids with Vorticity Diffusion. IEEE Transactions on Visualization and Computer Graphics, 2016.
diff --git a/SPlisHSPlasH/BoundaryModel.cpp b/SPlisHSPlasH/BoundaryModel.cpp
index 54f9352..1c07734 100644
--- a/SPlisHSPlasH/BoundaryModel.cpp
+++ b/SPlisHSPlasH/BoundaryModel.cpp
@@ -15,7 +15,6 @@ BoundaryModel::BoundaryModel() :
m_x(),
m_v(),
m_V(),
- m_boundaryPsi(),
m_forcePerThread(),
m_torquePerThread()
{
@@ -29,7 +28,6 @@ BoundaryModel::~BoundaryModel(void)
m_x.clear();
m_v.clear();
m_V.clear();
- m_boundaryPsi.clear();
m_forcePerThread.clear();
m_torquePerThread.clear();
@@ -58,7 +56,7 @@ void BoundaryModel::reset()
}
}
-void BoundaryModel::computeBoundaryPsi(const Real density0)
+void BoundaryModel::computeBoundaryVolume()
{
Simulation *sim = Simulation::getCurrent();
const unsigned int nFluids = sim->numberOfFluidModels();
@@ -83,7 +81,6 @@ void BoundaryModel::computeBoundaryPsi(const Real density0)
}
const Real volume = static_cast<Real>(1.0) / delta;
m_V[i] = volume;
- m_boundaryPsi[i] = density0 * volume;
}
}
}
@@ -94,7 +91,6 @@ void BoundaryModel::initModel(RigidBodyObject *rbo, const unsigned int numBounda
m_x.resize(numBoundaryParticles);
m_v.resize(numBoundaryParticles);
m_V.resize(numBoundaryParticles);
- m_boundaryPsi.resize(numBoundaryParticles);
#ifdef _OPENMP
const int maxThreads = omp_get_max_threads();
@@ -114,7 +110,6 @@ void BoundaryModel::initModel(RigidBodyObject *rbo, const unsigned int numBounda
m_x[i] = boundaryParticles[i];
m_v[i].setZero();
m_V[i] = 0.0;
- m_boundaryPsi[i] = 0.0;
}
}
m_rigidBody = rbo;
@@ -137,7 +132,6 @@ void BoundaryModel::performNeighborhoodSearchSort()
d.sort_field(&m_x[0]);
d.sort_field(&m_v[0]);
d.sort_field(&m_V[0]);
- d.sort_field(&m_boundaryPsi[0]);
m_sorted = true;
}
diff --git a/SPlisHSPlasH/BoundaryModel.h b/SPlisHSPlasH/BoundaryModel.h
index fbf3fa5..605c858 100644
--- a/SPlisHSPlasH/BoundaryModel.h
+++ b/SPlisHSPlasH/BoundaryModel.h
@@ -25,7 +25,6 @@ namespace SPH
std::vector<Vector3r> m_x;
std::vector<Vector3r> m_v;
std::vector<Real> m_V;
- std::vector<Real> m_boundaryPsi;
std::vector<Vector3r> m_forcePerThread;
std::vector<Vector3r> m_torquePerThread;
bool m_sorted;
@@ -34,7 +33,7 @@ namespace SPH
public:
unsigned int numberOfParticles() const { return static_cast<unsigned int>(m_x.size()); }
- void computeBoundaryPsi(const Real density0);
+ void computeBoundaryVolume();
virtual void reset();
@@ -119,21 +118,6 @@ namespace SPH
{
m_V[i] = val;
}
-
- FORCE_INLINE const Real& getBoundaryPsi(const unsigned int i) const
- {
- return m_boundaryPsi[i];
- }
-
- FORCE_INLINE Real& getBoundaryPsi(const unsigned int i)
- {
- return m_boundaryPsi[i];
- }
-
- FORCE_INLINE void setBoundaryPsi(const unsigned int i, const Real &val)
- {
- m_boundaryPsi[i] = val;
- }
};
}
diff --git a/SPlisHSPlasH/CMakeLists.txt b/SPlisHSPlasH/CMakeLists.txt
index 5f22c96..811066a 100644
--- a/SPlisHSPlasH/CMakeLists.txt
+++ b/SPlisHSPlasH/CMakeLists.txt
@@ -120,13 +120,27 @@ set(DRAG_SOURCE_FILES
Drag/DragForce_Macklin2014.cpp
)
+set(ELASTICITY_HEADER_FILES
+ Elasticity/ElasticityBase.h
+ Elasticity/Elasticity_Becker2009.h
+ Elasticity/Elasticity_Peer2018.h
+ )
+
+set(ELASTICITY_SOURCE_FILES
+ Elasticity/ElasticityBase.cpp
+ Elasticity/Elasticity_Becker2009.cpp
+ Elasticity/Elasticity_Peer2018.cpp
+ )
+
set(UTILS_HEADER_FILES
- Utilities/PoissonDiskSampling.h
+ Utilities/MathFunctions.h
Utilities/MatrixFreeSolver.h
+ Utilities/PoissonDiskSampling.h
Utilities/SceneLoader.h
)
set(UTILS_SOURCE_FILES
+ Utilities/MathFunctions.cpp
Utilities/PoissonDiskSampling.cpp
Utilities/SceneLoader.cpp
)
@@ -200,6 +214,9 @@ add_library(SPlisHSPlasH
${DRAG_HEADER_FILES}
${DRAG_SOURCE_FILES}
+ ${ELASTICITY_HEADER_FILES}
+ ${ELASTICITY_SOURCE_FILES}
+
${UTILS_HEADER_FILES}
${UTILS_SOURCE_FILES}
)
@@ -227,6 +244,8 @@ source_group("Header Files\\Vorticity" FILES ${VORTICITY_HEADER_FILES})
source_group("Source Files\\Vorticity" FILES ${VORTICITY_SOURCE_FILES})
source_group("Header Files\\Drag" FILES ${DRAG_HEADER_FILES})
source_group("Source Files\\Drag" FILES ${DRAG_SOURCE_FILES})
+source_group("Header Files\\Elasticity" FILES ${ELASTICITY_HEADER_FILES})
+source_group("Source Files\\Elasticity" FILES ${ELASTICITY_SOURCE_FILES})
source_group("Header Files\\Utils" FILES ${UTILS_HEADER_FILES})
source_group("Source Files\\Utils" FILES ${UTILS_SOURCE_FILES})
diff --git a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
index 28558f8..0b98dd4 100644
--- a/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
+++ b/SPlisHSPlasH/DFSPH/TimeStepDFSPH.cpp
@@ -28,10 +28,31 @@ TimeStepDFSPH::TimeStepDFSPH() :
m_enableDivergenceSolver = true;
m_maxIterationsV = 100;
m_maxErrorV = 0.1;
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "factor", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getFactor(fluidModelIndex, i); } });
+ model->addField({ "advected density", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDensityAdv(fluidModelIndex, i); } });
+ model->addField({ "kappa", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getKappa(fluidModelIndex, i); } });
+ model->addField({ "kappa_v", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getKappaV(fluidModelIndex, i); } });
+ }
}
TimeStepDFSPH::~TimeStepDFSPH(void)
{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("factor");
+ model->removeFieldByName("advected density");
+ model->removeFieldByName("kappa");
+ model->removeFieldByName("kappa_v");
+ }
}
void TimeStepDFSPH::initParameters()
@@ -206,7 +227,7 @@ void TimeStepDFSPH::warmstartPressureSolve(const unsigned int fluidModelIndex)
const Real invH2 = 1.0 / h2;
Simulation *sim = Simulation::getCurrent();
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const int numParticles = (int)model->numActiveParticles();
const unsigned int nFluids = sim->numberOfFluidModels();
@@ -290,7 +311,7 @@ void TimeStepDFSPH::pressureSolve()
for (unsigned int fluidModelIndex = 0; fluidModelIndex < nFluids; fluidModelIndex++)
{
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const int numParticles = (int)model->numActiveParticles();
#pragma omp parallel default(shared)
{
@@ -322,7 +343,7 @@ void TimeStepDFSPH::pressureSolve()
for (unsigned int i = 0; i < nFluids; i++)
{
FluidModel *model = sim->getFluidModel(i);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
avg_density_err = 0.0;
pressureSolveIteration(i, avg_density_err);
@@ -358,7 +379,7 @@ void TimeStepDFSPH::pressureSolveIteration(const unsigned int fluidModelIndex, R
{
Simulation *sim = Simulation::getCurrent();
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const int numParticles = (int)model->numActiveParticles();
const unsigned int nFluids = sim->numberOfFluidModels();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
@@ -441,7 +462,7 @@ void TimeStepDFSPH::warmstartDivergenceSolve(const unsigned int fluidModelIndex)
const Real invH = 1.0 / h;
Simulation *sim = Simulation::getCurrent();
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const int numParticles = (int)model->numActiveParticles();
const unsigned int nFluids = sim->numberOfFluidModels();
@@ -560,7 +581,7 @@ void TimeStepDFSPH::divergenceSolve()
for (unsigned int i = 0; i < nFluids; i++)
{
FluidModel *model = sim->getFluidModel(i);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
avg_density_err = 0.0;
divergenceSolveIteration(i, avg_density_err);
@@ -602,7 +623,7 @@ void TimeStepDFSPH::divergenceSolveIteration(const unsigned int fluidModelIndex,
{
Simulation *sim = Simulation::getCurrent();
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const int numParticles = (int)model->numActiveParticles();
const unsigned int nFluids = sim->numberOfFluidModels();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
diff --git a/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp b/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp
index 37374a9..80947b1 100644
--- a/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp
+++ b/SPlisHSPlasH/Drag/DragForce_Gissler2017.cpp
@@ -27,7 +27,7 @@ void DragForce_Gissler2017::step()
const Real radius = sim->getValue<Real>(Simulation::PARTICLE_RADIUS);
const Real diam = static_cast<Real>(2.0)*radius;
static const Real pi = static_cast<Real>(M_PI);
- const Real rho_l = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real rho_l = m_model->getDensity0();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
diff --git a/SPlisHSPlasH/Drag/DragForce_Macklin2014.cpp b/SPlisHSPlasH/Drag/DragForce_Macklin2014.cpp
index 8f0eda2..77b72e5 100644
--- a/SPlisHSPlasH/Drag/DragForce_Macklin2014.cpp
+++ b/SPlisHSPlasH/Drag/DragForce_Macklin2014.cpp
@@ -14,7 +14,7 @@ DragForce_Macklin2014::~DragForce_Macklin2014(void)
void DragForce_Macklin2014::step()
{
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const unsigned int numParticles = m_model->numActiveParticles();
diff --git a/SPlisHSPlasH/Elasticity/ElasticityBase.cpp b/SPlisHSPlasH/Elasticity/ElasticityBase.cpp
new file mode 100644
index 0000000..7426afb
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/ElasticityBase.cpp
@@ -0,0 +1,40 @@
+#include "ElasticityBase.h"
+
+using namespace SPH;
+using namespace GenParam;
+
+int ElasticityBase::YOUNGS_MODULUS = -1;
+int ElasticityBase::POISSON_RATIO = -1;
+
+
+ElasticityBase::ElasticityBase(FluidModel *model) :
+ NonPressureForceBase(model),
+ m_youngsModulus(100000.0),
+ m_poissonRatio(0.3)
+{
+}
+
+ElasticityBase::~ElasticityBase(void)
+{
+}
+
+
+void ElasticityBase::initParameters()
+{
+ NonPressureForceBase::initParameters();
+
+ YOUNGS_MODULUS = createNumericParameter("youngsModulus", "Young`s modulus", &m_youngsModulus);
+ setGroup(YOUNGS_MODULUS, "Elasticity");
+ setDescription(YOUNGS_MODULUS, "Stiffness of the elastic material");
+ RealParameter* rparam = static_cast<RealParameter*>(getParameter(YOUNGS_MODULUS));
+ rparam->setMinValue(0.0);
+
+ POISSON_RATIO = createNumericParameter("poissonsRatio", "Poisson`s ratio", &m_poissonRatio);
+ setGroup(POISSON_RATIO, "Elasticity");
+ setDescription(POISSON_RATIO, "Ratio of transversal expansion and axial compression");
+ rparam = static_cast<RealParameter*>(getParameter(POISSON_RATIO));
+ rparam->setMinValue(-1.0 + 1e-4);
+ rparam->setMaxValue(0.5 - 1e-4);
+}
+
+
diff --git a/SPlisHSPlasH/Elasticity/ElasticityBase.h b/SPlisHSPlasH/Elasticity/ElasticityBase.h
new file mode 100644
index 0000000..e1eae94
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/ElasticityBase.h
@@ -0,0 +1,29 @@
+#ifndef __ElasticityBase_h__
+#define __ElasticityBase_h__
+
+#include "SPlisHSPlasH/Common.h"
+#include "SPlisHSPlasH/FluidModel.h"
+#include "SPlisHSPlasH/NonPressureForceBase.h"
+
+namespace SPH
+{
+ /** \brief Base class for all elasticity methods.
+ */
+ class ElasticityBase : public NonPressureForceBase
+ {
+ protected:
+ Real m_youngsModulus;
+ Real m_poissonRatio;
+
+ virtual void initParameters();
+
+ public:
+ static int YOUNGS_MODULUS;
+ static int POISSON_RATIO;
+
+ ElasticityBase(FluidModel *model);
+ virtual ~ElasticityBase(void);
+ };
+}
+
+#endif
diff --git a/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.cpp b/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.cpp
new file mode 100644
index 0000000..a82b07f
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.cpp
@@ -0,0 +1,335 @@
+#include "Elasticity_Becker2009.h"
+#include "SPlisHSPlasH/Simulation.h"
+#include "SPlisHSPlasH/Utilities/MathFunctions.h"
+
+using namespace SPH;
+using namespace GenParam;
+
+int Elasticity_Becker2009::ALPHA = -1;
+
+
+Elasticity_Becker2009::Elasticity_Becker2009(FluidModel *model) :
+ ElasticityBase(model)
+{
+ const unsigned int numParticles = model->numActiveParticles();
+ m_restVolumes.resize(numParticles);
+ m_current_to_initial_index.resize(numParticles);
+ m_initial_to_current_index.resize(numParticles);
+ m_initialNeighbors.resize(numParticles);
+ m_rotations.resize(numParticles, Matrix3r::Identity());
+ m_stress.resize(numParticles);
+ m_F.resize(numParticles);
+ m_alpha = 0.0;
+
+ initValues();
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "rest volume", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &m_restVolumes[i]; } });
+ model->addField({ "rotation", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_rotations[i](0,0); } });
+ model->addField({ "stress", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_stress[i][0]; } });
+ model->addField({ "deformation gradient", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_F[i](0,0); } });
+ }
+}
+
+Elasticity_Becker2009::~Elasticity_Becker2009(void)
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("rest volume");
+ model->removeFieldByName("rotation");
+ model->removeFieldByName("stress");
+ model->removeFieldByName("deformation gradient");
+ }
+}
+
+void Elasticity_Becker2009::initParameters()
+{
+ ElasticityBase::initParameters();
+
+ ALPHA = createNumericParameter("alpha", "Zero-energy modes suppression", &m_alpha);
+ setGroup(ALPHA, "Elasticity");
+ setDescription(ALPHA, "Coefficent for zero-energy modes suppression method");
+ RealParameter *rparam = static_cast<RealParameter*>(getParameter(ALPHA));
+ rparam->setMinValue(0.0);
+}
+
+void Elasticity_Becker2009::initValues()
+{
+ Simulation *sim = Simulation::getCurrent();
+ sim->getNeighborhoodSearch()->find_neighbors();
+
+ FluidModel *model = m_model;
+ const unsigned int numParticles = model->numActiveParticles();
+ const unsigned int fluidModelIndex = model->getPointSetIndex();
+
+ // Store the neighbors in the reference configurations and
+ // compute the volume of each particle in rest state
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ m_current_to_initial_index[i] = i;
+ m_initial_to_current_index[i] = i;
+
+ // only neighbors in same phase will influence elasticity
+ const unsigned int numNeighbors = sim->numberOfNeighbors(fluidModelIndex, fluidModelIndex, i);
+ m_initialNeighbors[i].resize(numNeighbors);
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ m_initialNeighbors[i][j] = sim->getNeighbor(fluidModelIndex, fluidModelIndex, i, j);
+
+ // compute volume
+ Real density = model->getMass(i) * sim->W_zero();
+ const Vector3r &xi = model->getPosition(i);
+ forall_fluid_neighbors_in_same_phase(
+ density += model->getMass(neighborIndex) * sim->W(xi - xj);
+ )
+ m_restVolumes[i] = model->getMass(i) / density;
+ }
+ }
+}
+
+void Elasticity_Becker2009::step()
+{
+ computeRotations();
+ computeStress();
+ computeForces();
+}
+
+
+void Elasticity_Becker2009::reset()
+{
+ initValues();
+}
+
+void Elasticity_Becker2009::performNeighborhoodSearchSort()
+{
+ const unsigned int numPart = m_model->numActiveParticles();
+ if (numPart == 0)
+ return;
+
+ Simulation *sim = Simulation::getCurrent();
+ auto const& d = sim->getNeighborhoodSearch()->point_set(m_model->getPointSetIndex());
+ d.sort_field(&m_restVolumes[0]);
+ d.sort_field(&m_current_to_initial_index[0]);
+
+ for (unsigned int i = 0; i < numPart; i++)
+ m_initial_to_current_index[m_current_to_initial_index[i]] = i;
+}
+
+void Elasticity_Becker2009::computeRotations()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi = m_model->getPosition(i);
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+ Matrix3r Apq;
+ Apq.setZero();
+
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+ const Vector3r xj_xi = xj - xi;
+ const Vector3r xj_xi_0 = xj0 - xi0;
+ Apq += m_model->getMass(neighborIndex) * sim->W(xj_xi_0) * (xj_xi * xj_xi_0.transpose());
+ }
+
+// Vector3r sigma;
+// Matrix3r U, VT;
+// MathFunctions::svdWithInversionHandling(Apq, sigma, U, VT);
+// m_rotations[i] = U * VT;
+ Quaternionr q(m_rotations[i]);
+ MathFunctions::extractRotation(Apq, q, 10);
+ m_rotations[i] = q.matrix();
+ }
+ }
+}
+
+void Elasticity_Becker2009::computeStress()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+
+ // Elasticity tensor
+ Matrix6r C;
+ C.setZero();
+ const Real factor = m_youngsModulus / ((1.0 + m_poissonRatio)*(1.0 - 2.0 * m_poissonRatio));
+ C(0, 0) = C(1, 1) = C(2, 2) = factor * (1.0 - m_poissonRatio);
+ C(0, 1) = C(0, 2) = C(1, 0) = C(1, 2) = C(2, 0) = C(2, 1) = factor * (m_poissonRatio);
+ C(3, 3) = C(4, 4) = C(5, 5) = factor * 0.5*(1.0 - 2.0 * m_poissonRatio);
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi = m_model->getPosition(i);
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+
+ Matrix3r nablaU;
+ nablaU.setZero();
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+
+ const Vector3r xj_xi = xj - xi;
+ const Vector3r xj_xi_0 = xj0 - xi0;
+
+ const Vector3r uji = m_rotations[i].transpose() * xj_xi - xj_xi_0;
+ // subtract because kernel gradient is taken in direction of xji0 instead of xij0
+ nablaU -= (m_restVolumes[neighborIndex] * uji) * sim->gradW(xj_xi_0).transpose();
+ }
+ m_F[i] = nablaU + Matrix3r::Identity();
+
+ // compute Cauchy strain: epsilon = 0.5 (nabla u + nabla u^T)
+ Vector6r strain;
+ strain[0] = nablaU(0, 0); // \epsilon_{00}
+ strain[1] = nablaU(1, 1); // \epsilon_{11}
+ strain[2] = nablaU(2, 2); // \epsilon_{22}
+ strain[3] = 0.5 * (nablaU(0, 1) + nablaU(1, 0)); // \epsilon_{01}
+ strain[4] = 0.5 * (nablaU(0, 2) + nablaU(2, 0)); // \epsilon_{02}
+ strain[5] = 0.5 * (nablaU(1, 2) + nablaU(2, 1)); // \epsilon_{12}
+
+ // stress = C * epsilon
+ m_stress[i] = C * strain;
+ }
+ }
+}
+
+void Elasticity_Becker2009::computeForces()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+ Vector3r fi;
+ fi.setZero();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+
+ const Vector3r xj_xi_0 = xj0 - xi0;
+ const Vector3r gradW0 = sim->gradW(xj_xi_0);
+
+ const Vector3r dji = m_restVolumes[i] * gradW0;
+ const Vector3r dij = -m_restVolumes[neighborIndex] * gradW0;
+
+ Vector3r sdji, sdij;
+ symMatTimesVec(m_stress[neighborIndex], dji, sdji);
+ symMatTimesVec(m_stress[i], dij, sdij);
+
+ const Vector3r fij = -m_restVolumes[neighborIndex] * sdji;
+ const Vector3r fji = -m_restVolumes[i] * sdij;
+
+ fi += m_rotations[neighborIndex] * fij - m_rotations[i] * fji;
+ }
+ fi = 0.5*fi;
+
+ if (m_alpha != 0.0)
+ {
+ //////////////////////////////////////////////////////////////////////////
+ // Ganzenm�ller, G.C. 2015. An hourglass control algorithm for Lagrangian
+ // Smooth Particle Hydrodynamics. Computer Methods in Applied Mechanics and
+ // Engineering 286, 87�106.
+ //////////////////////////////////////////////////////////////////////////
+ Vector3r fi_hg;
+ fi_hg.setZero();
+ const Vector3r &xi = m_model->getPosition(i);
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+
+ // Note: Ganzenm�ller defines xij = xj-xi
+ const Vector3r xi_xj = -(xi - xj);
+ const Real xixj_l = xi_xj.norm();
+ if (xixj_l > 1.0e-6)
+ {
+ // Note: Ganzenm�ller defines xij = xj-xi
+ const Vector3r xi_xj_0 = -(xi0 - xj0);
+ const Real xixj0_l2 = xi_xj_0.squaredNorm();
+ const Real W0 = sim->W(xi_xj_0);
+
+ const Vector3r xij_i = m_F[i] * m_rotations[i] * xi_xj_0;
+ const Vector3r xji_j = -m_F[neighborIndex] * m_rotations[neighborIndex] * xi_xj_0;
+ const Vector3r epsilon_ij_i = xij_i - xi_xj;
+ const Vector3r epsilon_ji_j = xji_j + xi_xj;
+
+ const Real delta_ij_i = epsilon_ij_i.dot(xi_xj) / xixj_l;
+ const Real delta_ji_j = -epsilon_ji_j.dot(xi_xj) / xixj_l;
+
+ fi_hg -= m_restVolumes[neighborIndex] * W0 / xixj0_l2 * (delta_ij_i + delta_ji_j) * xi_xj / xixj_l;
+ }
+ }
+ fi_hg *= m_alpha * m_youngsModulus * m_restVolumes[i];
+ model->getAcceleration(i) += fi_hg / model->getMass(i);
+ }
+
+ // elastic acceleration
+ Vector3r &ai = model->getAcceleration(i);
+ ai += fi / model->getMass(i);
+ }
+ }
+}
+
diff --git a/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.h b/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.h
new file mode 100644
index 0000000..0d9ec4c
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/Elasticity_Becker2009.h
@@ -0,0 +1,59 @@
+#ifndef __Elasticity_Becker2009_h__
+#define __Elasticity_Becker2009_h__
+
+#include "SPlisHSPlasH/Common.h"
+#include "SPlisHSPlasH/FluidModel.h"
+#include "ElasticityBase.h"
+
+namespace SPH
+{
+ /** \brief This class implements the corotated SPH method for deformable solids introduced
+ * by Becker et al. \cite Becker:2009.
+ */
+ class Elasticity_Becker2009 : public ElasticityBase
+ {
+ protected:
+ // initial particle indices, used to access their original positions
+ std::vector<unsigned int> m_current_to_initial_index;
+ std::vector<unsigned int> m_initial_to_current_index;
+ // initial particle neighborhood
+ std::vector<std::vector<unsigned int>> m_initialNeighbors;
+ // volumes in rest configuration
+ std::vector<Real> m_restVolumes;
+ std::vector<Matrix3r> m_rotations;
+ std::vector<Vector6r> m_stress;
+ std::vector<Matrix3r> m_F;
+ Real m_alpha;
+
+ void initValues();
+ void computeRotations();
+ void computeStress();
+ void computeForces();
+
+ virtual void initParameters();
+
+ //////////////////////////////////////////////////////////////////////////
+ // multiplication of symmetric matrix, represented by a 6D vector, and a
+ // 3D vector
+ //////////////////////////////////////////////////////////////////////////
+ FORCE_INLINE void symMatTimesVec(const Vector6r & M, const Vector3r & v, Vector3r &res)
+ {
+ res[0] = M[0] * v[0] + M[3] * v[1] + M[4] * v[2];
+ res[1] = M[3] * v[0] + M[1] * v[1] + M[5] * v[2];
+ res[2] = M[4] * v[0] + M[5] * v[1] + M[2] * v[2];
+ }
+
+
+ public:
+ static int ALPHA;
+
+ Elasticity_Becker2009(FluidModel *model);
+ virtual ~Elasticity_Becker2009(void);
+
+ virtual void step();
+ virtual void reset();
+ virtual void performNeighborhoodSearchSort();
+ };
+}
+
+#endif
diff --git a/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.cpp b/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.cpp
new file mode 100644
index 0000000..866a51c
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.cpp
@@ -0,0 +1,562 @@
+#include "Elasticity_Peer2018.h"
+#include "SPlisHSPlasH/Simulation.h"
+#include "SPlisHSPlasH/Utilities/MathFunctions.h"
+#include "SPlisHSPlasH/TimeManager.h"
+#include "Utilities/Timing.h"
+#include "Utilities/Counting.h"
+
+using namespace SPH;
+using namespace GenParam;
+
+int Elasticity_Peer2018::ITERATIONS = -1;
+int Elasticity_Peer2018::MAX_ITERATIONS = -1;
+int Elasticity_Peer2018::MAX_ERROR = -1;
+int Elasticity_Peer2018::ALPHA = -1;
+
+
+Elasticity_Peer2018::Elasticity_Peer2018(FluidModel *model) :
+ ElasticityBase(model)
+{
+ const unsigned int numParticles = model->numActiveParticles();
+ m_restVolumes.resize(numParticles);
+ m_current_to_initial_index.resize(numParticles);
+ m_initial_to_current_index.resize(numParticles);
+ m_initialNeighbors.resize(numParticles);
+ m_rotations.resize(numParticles, Matrix3r::Identity());
+ m_stress.resize(numParticles);
+ m_L.resize(numParticles);
+ m_F.resize(numParticles);
+
+ m_iterations = 0;
+ m_maxIter = 100;
+ m_maxError = 1.0e-4;
+ m_alpha = 0.0;
+
+ initValues();
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "rest volume", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &m_restVolumes[i]; } });
+ model->addField({ "rotation", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_rotations[i](0,0); } });
+ model->addField({ "stress", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_stress[i][0]; } });
+ model->addField({ "deformation gradient", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_F[i](0,0); } });
+ model->addField({ "correction matrix", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_L[i](0,0); } });
+ }
+}
+
+Elasticity_Peer2018::~Elasticity_Peer2018(void)
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("rest volume");
+ model->removeFieldByName("rotation");
+ model->removeFieldByName("stress");
+ model->removeFieldByName("deformation gradient");
+ model->removeFieldByName("correction matrix");
+ }
+}
+
+
+void Elasticity_Peer2018::initParameters()
+{
+ ElasticityBase::initParameters();
+
+ ITERATIONS = createNumericParameter("elasticityIterations", "Iterations", &m_iterations);
+ setGroup(ITERATIONS, "Elasticity");
+ setDescription(ITERATIONS, "Iterations required by the elasticity solver.");
+ getParameter(ITERATIONS)->setReadOnly(true);
+
+ MAX_ITERATIONS = createNumericParameter("elasticityMaxIter", "Max. iterations (elasticity)", &m_maxIter);
+ setGroup(MAX_ITERATIONS, "Elasticity");
+ setDescription(MAX_ITERATIONS, "Coefficient for the elasticity force computation");
+ static_cast<NumericParameter<unsigned int>*>(getParameter(MAX_ITERATIONS))->setMinValue(1);
+
+ MAX_ERROR = createNumericParameter("elasticityMaxError", "Max. elasticity error", &m_maxError);
+ setGroup(MAX_ERROR, "Elasticity");
+ setDescription(MAX_ERROR, "Coefficient for the elasticity force computation");
+ RealParameter * rparam = static_cast<RealParameter*>(getParameter(MAX_ERROR));
+ rparam->setMinValue(1e-7);
+
+ ALPHA = createNumericParameter("alpha", "Zero-energy modes suppression", &m_alpha);
+ setGroup(ALPHA, "Elasticity");
+ setDescription(ALPHA, "Coefficent for zero-energy modes suppression method");
+ rparam = static_cast<RealParameter*>(getParameter(ALPHA));
+ rparam->setMinValue(0.0);
+}
+
+
+void Elasticity_Peer2018::initValues()
+{
+ Simulation *sim = Simulation::getCurrent();
+ sim->getNeighborhoodSearch()->find_neighbors();
+
+ FluidModel *model = m_model;
+ const unsigned int numParticles = model->numActiveParticles();
+ const unsigned int fluidModelIndex = model->getPointSetIndex();
+
+ // Store the neighbors in the reference configurations and
+ // compute the volume of each particle in rest state
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ m_current_to_initial_index[i] = i;
+ m_initial_to_current_index[i] = i;
+
+ // only neighbors in same phase will influence elasticity
+ const unsigned int numNeighbors = sim->numberOfNeighbors(fluidModelIndex, fluidModelIndex, i);
+ m_initialNeighbors[i].resize(numNeighbors);
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ m_initialNeighbors[i][j] = sim->getNeighbor(fluidModelIndex, fluidModelIndex, i, j);
+
+ // compute volume
+ Real density = model->getMass(i) * sim->W_zero();
+ const Vector3r &xi = model->getPosition(i);
+ forall_fluid_neighbors_in_same_phase(
+ density += model->getMass(neighborIndex) * sim->W(xi - xj);
+ )
+ m_restVolumes[i] = model->getMass(i) / density;
+ m_rotations[i].setIdentity();
+ }
+ }
+
+ computeMatrixL();
+}
+
+
+void Elasticity_Peer2018::step()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const Real dt = TimeManager::getCurrent()->getTimeStepSize();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Init linear system solver and preconditioner
+ //////////////////////////////////////////////////////////////////////////
+ MatrixReplacement A(3 * m_model->numActiveParticles(), matrixVecProd, (void*)this);
+ //m_solver.preconditioner().init(m_model->numActiveParticles(), diagonalMatrixElement, (void*)this);
+
+ m_solver.setTolerance(m_maxError);
+ m_solver.setMaxIterations(m_maxIter);
+ m_solver.compute(A);
+
+ VectorXr b(3 * numParticles);
+ VectorXr x(3 * numParticles);
+ VectorXr g(3 * numParticles);
+
+ computeRotations();
+ computeRHS(b);
+
+ // warmstart
+ #pragma omp parallel for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ g.segment<3>(3 * i) = m_model->getVelocity(i) + dt * m_model->getAcceleration(i);
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ // Solve linear system
+ //////////////////////////////////////////////////////////////////////////
+ START_TIMING("Elasticity - CG solve");
+ x = m_solver.solveWithGuess(b, g);
+ m_iterations = (int)m_solver.iterations();
+ STOP_TIMING_AVG;
+ INCREASE_COUNTER("Elasticity - CG iterations", static_cast<Real>(m_iterations));
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ Vector3r &ai = m_model->getAcceleration(i);
+ ai += (1.0 / dt) * (x.segment<3>(3 * i) - m_model->getVelocity(i));
+ }
+ }
+}
+
+
+void Elasticity_Peer2018::reset()
+{
+ initValues();
+}
+
+void Elasticity_Peer2018::performNeighborhoodSearchSort()
+{
+ const unsigned int numPart = m_model->numActiveParticles();
+ if (numPart == 0)
+ return;
+
+ Simulation *sim = Simulation::getCurrent();
+ auto const& d = sim->getNeighborhoodSearch()->point_set(m_model->getPointSetIndex());
+ d.sort_field(&m_restVolumes[0]);
+ d.sort_field(&m_rotations[0]);
+ d.sort_field(&m_current_to_initial_index[0]);
+ d.sort_field(&m_L[0]);
+
+ for (unsigned int i = 0; i < numPart; i++)
+ m_initial_to_current_index[m_current_to_initial_index[i]] = i;
+}
+
+void Elasticity_Peer2018::computeRotations()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi = m_model->getPosition(i);
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+ Matrix3r F;
+ F.setZero();
+
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+ const Vector3r xj_xi = xj - xi;
+ const Vector3r xi_xj_0 = xi0 - xj0;
+ const Vector3r correctedKernel = m_L[i] * sim->gradW(xi_xj_0);
+ F += m_restVolumes[neighborIndex] * xj_xi * correctedKernel.transpose();
+ }
+
+// Vector3r sigma;
+// Matrix3r U, VT;
+// MathFunctions::svdWithInversionHandling(F, sigma, U, VT);
+// m_rotations[i] = U * VT;
+ Quaternionr q(m_rotations[i]);
+ MathFunctions::extractRotation(F, q, 10);
+ m_rotations[i] = q.matrix();
+ }
+ }
+}
+
+void Elasticity_Peer2018::computeMatrixL()
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+ Matrix3r L;
+ L.setZero();
+
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+ const Vector3r xj_xi_0 = xj0 - xi0;
+ const Vector3r gradW = sim->gradW(xj_xi_0);
+
+ // minus because gradW(xij0) == -gradW(xji0)
+ L -= m_restVolumes[neighborIndex] * gradW * xj_xi_0.transpose();
+ }
+
+ bool invertible = false;
+ L.computeInverseWithCheck(m_L[i], invertible, 1e-9);
+ if (!invertible)
+ {
+ //MathFunctions::pseudoInverse(L, m_L[i]);
+ m_L[i] = Matrix3r::Zero();
+ }
+ }
+ }
+}
+
+
+void Elasticity_Peer2018::computeRHS(VectorXr & rhs)
+{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int numParticles = m_model->numActiveParticles();
+ const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ FluidModel *model = m_model;
+ const Real dt = TimeManager::getCurrent()->getTimeStepSize();
+
+ Real mu = m_youngsModulus / (static_cast<Real>(2.0) * (static_cast<Real>(1.0) + m_poissonRatio));
+ Real lambda = m_youngsModulus * m_poissonRatio / ((static_cast<Real>(1.0) + m_poissonRatio) * (static_cast<Real>(1.0) - static_cast<Real>(2.0) * m_poissonRatio));
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi = m_model->getPosition(i);
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // compute corotated deformation gradient (Eq. 18)
+ //////////////////////////////////////////////////////////////////////////
+ m_F[i].setZero();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+ const Vector3r xj_xi = xj - xi;
+ const Vector3r xi_xj_0 = xi0 - xj0;
+ const Vector3r correctedRotatedKernel = m_rotations[i] * m_L[i] * sim->gradW(xi_xj_0);
+ m_F[i] += m_restVolumes[neighborIndex] * xj_xi * correctedRotatedKernel.transpose();
+ }
+
+ //////////////////////////////////////////////////////////////////////////
+ // compute Cauchy strain: epsilon = 0.5 (F + F^T) - I
+ //////////////////////////////////////////////////////////////////////////
+ Vector6r strain;
+ strain[0] = m_F[i](0, 0) - 1.0; // \epsilon_{00}
+ strain[1] = m_F[i](1, 1) - 1.0; // \epsilon_{11}
+ strain[2] = m_F[i](2, 2) - 1.0; // \epsilon_{22}
+ strain[3] = 0.5 * (m_F[i](0, 1) + m_F[i](1, 0)); // \epsilon_{01}
+ strain[4] = 0.5 * (m_F[i](0, 2) + m_F[i](2, 0)); // \epsilon_{02}
+ strain[5] = 0.5 * (m_F[i](1, 2) + m_F[i](2, 1)); // \epsilon_{12}
+
+ //////////////////////////////////////////////////////////////////////////
+ // First Piola Kirchhoff stress = 2 mu epsilon + lambda trace(epsilon) I
+ //////////////////////////////////////////////////////////////////////////
+ const Real trace = strain[0] + strain[1] + strain[2];
+ const Real ltrace = lambda*trace;
+ m_stress[i] = strain * 2.0*mu;
+ m_stress[i][0] += ltrace;
+ m_stress[i][1] += ltrace;
+ m_stress[i][2] += ltrace;
+ }
+ }
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = m_current_to_initial_index[i];
+ const Vector3r &xi0 = m_model->getPosition0(i0);
+ const size_t numNeighbors = m_initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Compute elastic force
+ //////////////////////////////////////////////////////////////////////////
+ Vector3r force;
+ force.setZero();
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+ const Vector3r xi_xj_0 = xi0 - xj0;
+ const Vector3r correctedRotatedKernel_i = m_rotations[i] * m_L[i] * sim->gradW(xi_xj_0);
+ const Vector3r correctedRotatedKernel_j = -m_rotations[neighborIndex] * m_L[neighborIndex] * sim->gradW(xi_xj_0);
+ Vector3r PWi, PWj;
+ symMatTimesVec(m_stress[i], correctedRotatedKernel_i, PWi);
+ symMatTimesVec(m_stress[neighborIndex], correctedRotatedKernel_j, PWj);
+ force += m_restVolumes[i] * m_restVolumes[neighborIndex] * (PWi - PWj);
+ }
+
+ if (m_alpha != 0.0)
+ {
+ //////////////////////////////////////////////////////////////////////////
+ // Ganzenm�ller, G.C. 2015. An hourglass control algorithm for Lagrangian
+ // Smooth Particle Hydrodynamics. Computer Methods in Applied Mechanics and
+ // Engineering 286, 87�106.
+ //////////////////////////////////////////////////////////////////////////
+ Vector3r fi_hg;
+ fi_hg.setZero();
+ const Vector3r &xi = m_model->getPosition(i);
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
+
+ const Vector3r &xj = model->getPosition(neighborIndex);
+ const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
+
+ // Note: Ganzenm�ller defines xij = xj-xi
+ const Vector3r xi_xj = -(xi - xj);
+ const Real xixj_l = xi_xj.norm();
+ if (xixj_l > 1.0e-6)
+ {
+ // Note: Ganzenm�ller defines xij = xj-xi
+ const Vector3r xi_xj_0 = -(xi0 - xj0);
+ const Real xixj0_l2 = xi_xj_0.squaredNorm();
+ const Real W0 = sim->W(xi_xj_0);
+
+ const Vector3r xij_i = m_F[i] * m_rotations[i] * xi_xj_0;
+ const Vector3r xji_j = -m_F[neighborIndex] * m_rotations[neighborIndex] * xi_xj_0;
+ const Vector3r epsilon_ij_i = xij_i - xi_xj;
+ const Vector3r epsilon_ji_j = xji_j + xi_xj;
+
+ const Real delta_ij_i = epsilon_ij_i.dot(xi_xj) / xixj_l;
+ const Real delta_ji_j = -epsilon_ji_j.dot(xi_xj) / xixj_l;
+
+ fi_hg -= m_restVolumes[neighborIndex] * W0 / xixj0_l2 * (delta_ij_i + delta_ji_j) * xi_xj / xixj_l;
+ }
+ }
+ fi_hg *= m_alpha * m_youngsModulus * m_restVolumes[i];
+ model->getAcceleration(i) += fi_hg / model->getMass(i);
+ }
+
+ rhs.segment<3>(3 * i) = model->getVelocity(i) + dt * (model->getAcceleration(i) + 1.0 / model->getMass(i) * force);
+ }
+ }
+}
+
+void Elasticity_Peer2018::matrixVecProd(const Real* vec, Real *result, void *userData)
+{
+ Simulation *sim = Simulation::getCurrent();
+ Elasticity_Peer2018 * elasticity = static_cast<Elasticity_Peer2018*>(userData);
+ FluidModel *model = elasticity->getModel();
+ const unsigned int numParticles = model->numActiveParticles();
+ const Real dt = TimeManager::getCurrent()->getTimeStepSize();
+
+ const auto ¤t_to_initial_index = elasticity->m_current_to_initial_index;
+ const auto &initial_to_current_index = elasticity->m_initial_to_current_index;
+ const auto &initialNeighbors = elasticity->m_initialNeighbors;
+ const auto &restVolumes = elasticity->m_restVolumes;
+ const auto &rotations = elasticity->m_rotations;
+ const auto &L = elasticity->m_L;
+ auto &stress = elasticity->m_stress;
+ const Real youngsModulus = elasticity->m_youngsModulus;
+ const Real poissonRatio = elasticity->m_poissonRatio;
+
+ Real mu = youngsModulus / (static_cast<Real>(2.0) * (static_cast<Real>(1.0) + poissonRatio));
+ Real lambda = youngsModulus * poissonRatio / ((static_cast<Real>(1.0) + poissonRatio) * (static_cast<Real>(1.0) - static_cast<Real>(2.0) * poissonRatio));
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = current_to_initial_index[i];
+ const Vector3r &pi = Eigen::Map<const Vector3r>(&vec[3 * i], 3);
+ const Vector3r &xi0 = model->getPosition0(i0);
+ const size_t numNeighbors = initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // compute corotated deformation gradient (Eq. 18)
+ //////////////////////////////////////////////////////////////////////////
+ Matrix3r nablaU;
+ nablaU.setZero();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Fluid
+ //////////////////////////////////////////////////////////////////////////
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = initial_to_current_index[initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = initialNeighbors[i0][j];
+
+ const Vector3r &pj = Eigen::Map<const Vector3r>(&vec[3 * neighborIndex], 3);
+ const Vector3r &xj0 = model->getPosition0(neighborIndex0);
+ const Vector3r pj_pi = pj - pi;
+ const Vector3r xi_xj_0 = xi0 - xj0;
+ const Vector3r correctedRotatedKernel = rotations[i] * L[i] * sim->gradW(xi_xj_0);
+ nablaU += restVolumes[neighborIndex] * pj_pi * correctedRotatedKernel.transpose();
+ }
+ nablaU *= dt;
+
+ //////////////////////////////////////////////////////////////////////////
+ // compute Cauchy strain: epsilon = 0.5 (nablaU + nablaU^T)
+ //////////////////////////////////////////////////////////////////////////
+ Vector6r strain;
+ strain[0] = nablaU(0, 0); // \epsilon_{00}
+ strain[1] = nablaU(1, 1); // \epsilon_{11}
+ strain[2] = nablaU(2, 2); // \epsilon_{22}
+ strain[3] = 0.5 * (nablaU(0, 1) + nablaU(1, 0)); // \epsilon_{01}
+ strain[4] = 0.5 * (nablaU(0, 2) + nablaU(2, 0)); // \epsilon_{02}
+ strain[5] = 0.5 * (nablaU(1, 2) + nablaU(2, 1)); // \epsilon_{12}
+
+ //////////////////////////////////////////////////////////////////////////
+ // First Piola Kirchhoff stress = 2 mu epsilon + lambda trace(epsilon) I
+ //////////////////////////////////////////////////////////////////////////
+ const Real trace = strain[0] + strain[1] + strain[2];
+ const Real ltrace = lambda*trace;
+ stress[i] = strain * 2.0*mu;
+ stress[i][0] += ltrace;
+ stress[i][1] += ltrace;
+ stress[i][2] += ltrace;
+ }
+ }
+
+ #pragma omp parallel default(shared)
+ {
+ #pragma omp for schedule(static)
+ for (int i = 0; i < (int)numParticles; i++)
+ {
+ const unsigned int i0 = current_to_initial_index[i];
+ const Vector3r &xi0 = model->getPosition0(i0);
+ const size_t numNeighbors = initialNeighbors[i0].size();
+
+ //////////////////////////////////////////////////////////////////////////
+ // Compute elastic force
+ //////////////////////////////////////////////////////////////////////////
+ Vector3r force;
+ force.setZero();
+ for (unsigned int j = 0; j < numNeighbors; j++)
+ {
+ const unsigned int neighborIndex = initial_to_current_index[initialNeighbors[i0][j]];
+ // get initial neighbor index considering the current particle order
+ const unsigned int neighborIndex0 = initialNeighbors[i0][j];
+
+ const Vector3r &xj0 = model->getPosition0(neighborIndex0);
+ const Vector3r xi_xj_0 = xi0 - xj0;
+ const Vector3r correctedRotatedKernel_i = rotations[i] * L[i] * sim->gradW(xi_xj_0);
+ const Vector3r correctedRotatedKernel_j = -rotations[neighborIndex] * L[neighborIndex] * sim->gradW(xi_xj_0);
+ Vector3r PWi, PWj;
+ elasticity->symMatTimesVec(stress[i], correctedRotatedKernel_i, PWi);
+ elasticity->symMatTimesVec(stress[neighborIndex], correctedRotatedKernel_j, PWj);
+ force += restVolumes[i] * restVolumes[neighborIndex] * (PWi - PWj);
+ }
+
+ const Real factor = dt / model->getMass(i);
+ result[3 * i] = vec[3 * i] - factor * force[0];
+ result[3 * i + 1] = vec[3 * i + 1] - factor * force[1];
+ result[3 * i + 2] = vec[3 * i + 2] - factor * force[2];
+ }
+ }
+}
diff --git a/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h b/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h
new file mode 100644
index 0000000..b167df4
--- /dev/null
+++ b/SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h
@@ -0,0 +1,73 @@
+#ifndef __Elasticity_Peer2018_h__
+#define __Elasticity_Peer2018_h__
+
+#include "SPlisHSPlasH/Common.h"
+#include "SPlisHSPlasH/FluidModel.h"
+#include "ElasticityBase.h"
+#include "SPlisHSPlasH/Utilities/MatrixFreeSolver.h"
+
+namespace SPH
+{
+ /** \brief This class implements the implicit SPH formulation for
+ * incompressible linearly elastic solids introduced
+ * by Peer et al. \cite PGBT17.
+ */
+ class Elasticity_Peer2018 : public ElasticityBase
+ {
+ protected:
+ typedef Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> Solver;
+
+ // initial particle indices, used to access their original positions
+ std::vector<unsigned int> m_current_to_initial_index;
+ std::vector<unsigned int> m_initial_to_current_index;
+ // initial particle neighborhood
+ std::vector<std::vector<unsigned int>> m_initialNeighbors;
+ // volumes in rest configuration
+ std::vector<Real> m_restVolumes;
+ std::vector<Matrix3r> m_rotations;
+ std::vector<Vector6r> m_stress;
+ std::vector<Matrix3r> m_L;
+ std::vector<Matrix3r> m_F;
+ unsigned int m_iterations;
+ unsigned int m_maxIter;
+ Real m_maxError;
+ Real m_alpha;
+ Solver m_solver;
+
+ void initValues();
+ void computeMatrixL();
+ void computeRotations();
+ void computeRHS(VectorXr & rhs);
+
+ virtual void initParameters();
+
+ //////////////////////////////////////////////////////////////////////////
+ // multiplication of symmetric matrix, represented by a 6D vector, and a
+ // 3D vector
+ //////////////////////////////////////////////////////////////////////////
+ FORCE_INLINE void symMatTimesVec(const Vector6r & M, const Vector3r & v, Vector3r &res)
+ {
+ res[0] = M[0] * v[0] + M[3] * v[1] + M[4] * v[2];
+ res[1] = M[3] * v[0] + M[1] * v[1] + M[5] * v[2];
+ res[2] = M[4] * v[0] + M[5] * v[1] + M[2] * v[2];
+ }
+
+
+ public:
+ static int ITERATIONS;
+ static int MAX_ITERATIONS;
+ static int MAX_ERROR;
+ static int ALPHA;
+
+ Elasticity_Peer2018(FluidModel *model);
+ virtual ~Elasticity_Peer2018(void);
+
+ virtual void step();
+ virtual void reset();
+ virtual void performNeighborhoodSearchSort();
+
+ static void matrixVecProd(const Real* vec, Real *result, void *userData);
+ };
+}
+
+#endif
diff --git a/SPlisHSPlasH/FluidModel.cpp b/SPlisHSPlasH/FluidModel.cpp
index 6aed189..8cad952 100644
--- a/SPlisHSPlasH/FluidModel.cpp
+++ b/SPlisHSPlasH/FluidModel.cpp
@@ -25,6 +25,9 @@
#include "Vorticity/MicropolarModel_Bender2017.h"
#include "Drag/DragForce_Gissler2017.h"
#include "Drag/DragForce_Macklin2014.h"
+#include "Elasticity/ElasticityBase.h"
+#include "Elasticity/Elasticity_Becker2009.h"
+#include "Elasticity/Elasticity_Peer2018.h"
using namespace SPH;
@@ -37,6 +40,7 @@ int FluidModel::DRAG_METHOD = -1;
int FluidModel::SURFACE_TENSION_METHOD = -1;
int FluidModel::VISCOSITY_METHOD = -1;
int FluidModel::VORTICITY_METHOD = -1;
+int FluidModel::ELASTICITY_METHOD = -1;
int FluidModel::ENUM_DRAG_NONE = -1;
int FluidModel::ENUM_DRAG_MACKLIN2014 = -1;
int FluidModel::ENUM_DRAG_GISSLER2017 = -1;
@@ -55,6 +59,9 @@ int FluidModel::ENUM_VISCOSITY_WEILER2018 = -1;
int FluidModel::ENUM_VORTICITY_NONE = -1;
int FluidModel::ENUM_VORTICITY_MICROPOLAR = -1;
int FluidModel::ENUM_VORTICITY_VC = -1;
+int FluidModel::ENUM_ELASTICITY_NONE = -1;
+int FluidModel::ENUM_ELASTICITY_BECKER2009 = -1;
+int FluidModel::ENUM_ELASTICITY_PEER2018 = -1;
FluidModel::FluidModel() :
@@ -82,15 +89,27 @@ FluidModel::FluidModel() :
m_surfaceTensionMethodChanged = nullptr;
m_viscosityMethodChanged = nullptr;
m_vorticityMethodChanged = nullptr;
+ m_elasticityMethod = ElasticityMethods::None;
+ m_elasticity = nullptr;
+ m_elasticityMethodChanged = nullptr;
+
+ addField({ "position", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &getPosition(i)[0]; } });
+ addField({ "velocity", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &getVelocity(i)[0]; } });
+ addField({ "density", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &getDensity(i); } });
}
FluidModel::~FluidModel(void)
{
+ removeFieldByName("position");
+ removeFieldByName("velocity");
+ removeFieldByName("density");
+
delete m_emitterSystem;
delete m_surfaceTension;
delete m_drag;
delete m_vorticity;
delete m_viscosity;
+ delete m_elasticity;
releaseFluidParticles();
}
@@ -171,6 +190,16 @@ void FluidModel::initParameters()
enumParam->addEnumValue("None", ENUM_VORTICITY_NONE);
enumParam->addEnumValue("Micropolar model", ENUM_VORTICITY_MICROPOLAR);
enumParam->addEnumValue("Vorticity confinement", ENUM_VORTICITY_VC);
+
+ ParameterBase::GetFunc<int> getElasticityFct = std::bind(&FluidModel::getElasticityMethod, this);
+ ParameterBase::SetFunc<int> setElasticityFct = std::bind(&FluidModel::setElasticityMethod, this, std::placeholders::_1);
+ ELASTICITY_METHOD = createEnumParameter("elasticityMethod", "Elasticity method", getElasticityFct, setElasticityFct);
+ setGroup(ELASTICITY_METHOD, "Elasticity");
+ setDescription(ELASTICITY_METHOD, "Method to compute elastic forces.");
+ enumParam = static_cast<EnumParameter*>(getParameter(ELASTICITY_METHOD));
+ enumParam->addEnumValue("None", ENUM_ELASTICITY_NONE);
+ enumParam->addEnumValue("Becker et al. 2009", ENUM_ELASTICITY_BECKER2009);
+ enumParam->addEnumValue("Peer et al. 2018", ENUM_ELASTICITY_PEER2018);
}
@@ -201,6 +230,8 @@ void FluidModel::reset()
m_vorticity->reset();
if (m_drag)
m_drag->reset();
+ if (m_elasticity)
+ m_elasticity->reset();
m_emitterSystem->reset();
}
@@ -305,15 +336,29 @@ void FluidModel::performNeighborhoodSearchSort()
m_vorticity->performNeighborhoodSearchSort();
if (m_drag)
m_drag->performNeighborhoodSearchSort();
+ if (m_elasticity)
+ m_elasticity->performNeighborhoodSearchSort();
}
void SPH::FluidModel::setDensity0(const Real v)
{
m_density0 = v;
initMasses();
- Simulation::getCurrent()->updateBoundaryPsi();
}
+const SPH::FieldDescription & SPH::FluidModel::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 FluidModel::setNumActiveParticles(const unsigned int num)
{
@@ -345,6 +390,11 @@ void FluidModel::setVorticityMethodChangedCallback(std::function<void()> const&
m_vorticityMethodChanged = callBackFct;
}
+void FluidModel::setElasticityMethodChangedCallback(std::function<void()> const& callBackFct)
+{
+ m_elasticityMethodChanged = callBackFct;
+}
+
void FluidModel::computeSurfaceTension()
{
if (m_surfaceTension)
@@ -369,6 +419,12 @@ void FluidModel::computeDragForce()
m_drag->step();
}
+void FluidModel::computeElasticity()
+{
+ if (m_elasticity)
+ m_elasticity->step();
+}
+
void FluidModel::emittedParticles(const unsigned int startIndex)
{
if (m_viscosity)
@@ -379,6 +435,8 @@ void FluidModel::emittedParticles(const unsigned int startIndex)
m_vorticity->emittedParticles(startIndex);
if (m_drag)
m_drag->emittedParticles(startIndex);
+ if (m_elasticity)
+ m_elasticity->emittedParticles(startIndex);
}
void FluidModel::setSurfaceTensionMethod(const int val)
@@ -495,3 +553,49 @@ void FluidModel::setDragMethod(const int val)
if (m_dragMethodChanged != nullptr)
m_dragMethodChanged();
}
+
+void FluidModel::setElasticityMethod(const int val)
+{
+ ElasticityMethods em = static_cast<ElasticityMethods>(val);
+ if ((em < ElasticityMethods::None) || (em >= ElasticityMethods::NumElasticityMethods))
+ em = ElasticityMethods::None;
+
+ if (em == m_elasticityMethod)
+ return;
+
+ delete m_elasticity;
+ m_elasticity = nullptr;
+
+ m_elasticityMethod = em;
+
+ if (m_elasticityMethod == ElasticityMethods::Becker2009)
+ m_elasticity = new Elasticity_Becker2009(this);
+ else if (m_elasticityMethod == ElasticityMethods::Peer2018)
+ m_elasticity = new Elasticity_Peer2018(this);
+
+ if (m_elasticity != nullptr)
+ m_elasticity->init();
+
+ if (m_elasticityMethodChanged != nullptr)
+ m_elasticityMethodChanged();
+}
+
+
+void FluidModel::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); });
+}
+
+void FluidModel::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;
+ }
+ }
+}
+
diff --git a/SPlisHSPlasH/FluidModel.h b/SPlisHSPlasH/FluidModel.h
index e9f1b05..b8588c2 100644
--- a/SPlisHSPlasH/FluidModel.h
+++ b/SPlisHSPlasH/FluidModel.h
@@ -15,12 +15,23 @@ namespace SPH
class SurfaceTensionBase;
class VorticityBase;
class DragBase;
+ class ElasticityBase;
class EmitterSystem;
+ enum FieldType { Scalar = 0, Vector3, Vector6, Matrix3, Matrix6 };
+ struct FieldDescription
+ {
+ std::string name;
+ FieldType type;
+ // getFct(particleIndex)
+ std::function<Real*(const unsigned int)> getFct;
+ };
+
enum class SurfaceTensionMethods { None = 0, Becker2007, Akinci2013, He2014, NumSurfaceTensionMethods };
enum class ViscosityMethods { None = 0, Standard, XSPH, Bender2017, Peer2015, Peer2016, Takahashi2015, Weiler2018, NumViscosityMethods };
enum class VorticityMethods { None = 0, Micropolar, VorticityConfinement, NumVorticityMethods };
enum class DragMethods { None = 0, Macklin2014, Gissler2017, NumDragMethods };
+ enum class ElasticityMethods { None = 0, Becker2009, Peer2018, NumElasticityMethods };
/** \brief The fluid model stores the particle and simulation information
*/
@@ -35,6 +46,7 @@ namespace SPH
static int SURFACE_TENSION_METHOD;
static int VISCOSITY_METHOD;
static int VORTICITY_METHOD;
+ static int ELASTICITY_METHOD;
static int ENUM_DRAG_NONE;
static int ENUM_DRAG_MACKLIN2014;
@@ -58,7 +70,10 @@ namespace SPH
static int ENUM_VORTICITY_MICROPOLAR;
static int ENUM_VORTICITY_VC;
-
+ static int ENUM_ELASTICITY_NONE;
+ static int ENUM_ELASTICITY_BECKER2009;
+ static int ENUM_ELASTICITY_PEER2018;
+
FluidModel();
virtual ~FluidModel();
@@ -89,11 +104,15 @@ namespace SPH
VorticityBase *m_vorticity;
DragMethods m_dragMethod;
DragBase *m_drag;
+ ElasticityMethods m_elasticityMethod;
+ ElasticityBase *m_elasticity;
+ std::vector<FieldDescription> m_fields;
std::function<void()> m_dragMethodChanged;
std::function<void()> m_surfaceTensionMethodChanged;
std::function<void()> m_viscosityMethodChanged;
std::function<void()> m_vorticityMethodChanged;
+ std::function<void()> m_elasticityMethodChanged;
Real m_density0;
unsigned int m_pointSetIndex;
@@ -120,6 +139,13 @@ namespace SPH
unsigned int getPointSetIndex() const { return m_pointSetIndex; }
+ void addField(const FieldDescription &field);
+ const std::vector<FieldDescription> &getFields() { return m_fields; }
+ const FieldDescription &getField(const unsigned int i) { return m_fields[i]; }
+ const FieldDescription &getField(const std::string &name);
+ const unsigned int numberOfFields() { return static_cast<unsigned int>(m_fields.size()); }
+ void removeFieldByName(const std::string &fieldName);
+
void setNumActiveParticles(const unsigned int num);
unsigned int numberOfParticles() const { return static_cast<unsigned int>(m_x.size()); }
@@ -147,21 +173,26 @@ namespace SPH
void setVorticityMethod(const int val);
int getDragMethod() const { return static_cast<int>(m_dragMethod); }
void setDragMethod(const int val);
+ int getElasticityMethod() const { return static_cast<int>(m_elasticityMethod); }
+ void setElasticityMethod(const int val);
SurfaceTensionBase *getSurfaceTensionBase() { return m_surfaceTension; }
ViscosityBase *getViscosityBase() { return m_viscosity; }
VorticityBase *getVorticityBase() { return m_vorticity; }
DragBase *getDragBase() { return m_drag; }
+ ElasticityBase *getElasticityBase() { return m_elasticity; }
void setDragMethodChangedCallback(std::function<void()> const& callBackFct);
void setSurfaceMethodChangedCallback(std::function<void()> const& callBackFct);
void setViscosityMethodChangedCallback(std::function<void()> const& callBackFct);
void setVorticityMethodChangedCallback(std::function<void()> const& callBackFct);
+ void setElasticityMethodChangedCallback(std::function<void()> const& callBackFct);
void computeSurfaceTension();
void computeViscosity();
void computeVorticity();
void computeDragForce();
+ void computeElasticity();
FORCE_INLINE Vector3r &getPosition0(const unsigned int i)
diff --git a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
index efae128..1c97210 100644
--- a/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
+++ b/SPlisHSPlasH/IISPH/TimeStepIISPH.cpp
@@ -14,10 +14,37 @@ TimeStepIISPH::TimeStepIISPH() :
{
m_simulationData.init();
m_counter = 0;
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "a_ii", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getAii(fluidModelIndex, i); } });
+ model->addField({ "d_ii", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDii(fluidModelIndex, i)[0]; } });
+ model->addField({ "d_ij*p_j", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDij_pj(fluidModelIndex, i)[0]; } });
+ model->addField({ "pressure", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressure(fluidModelIndex, i); } });
+ //model->addField({ "last pressure", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getLastPressure(fluidModelIndex, i); } });
+ model->addField({ "advected density", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDensityAdv(fluidModelIndex, i); } });
+ model->addField({ "pressure acceleration", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressureAccel(fluidModelIndex, i)[0]; } });
+ }
}
TimeStepIISPH::~TimeStepIISPH(void)
{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("a_ii");
+ model->removeFieldByName("d_ii");
+ model->removeFieldByName("d_ij*p_j");
+ model->removeFieldByName("pressure");
+ //model->removeFieldByName("last pressure");
+ model->removeFieldByName("advected density");
+ model->removeFieldByName("pressure acceleration");
+ }
}
void TimeStepIISPH::step()
@@ -187,7 +214,7 @@ void TimeStepIISPH::pressureSolve()
for (unsigned int i = 0; i < nFluids; i++)
{
FluidModel *model = sim->getFluidModel(i);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
avg_density_err = 0.0;
pressureSolveIteration(i, avg_density_err);
@@ -208,7 +235,7 @@ void TimeStepIISPH::pressureSolveIteration(const unsigned int fluidModelIndex, R
const unsigned int numParticles = model->numActiveParticles();
const unsigned int nFluids = sim->numberOfFluidModels();
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
const Real h2 = h*h;
const Real omega = 0.5;
diff --git a/SPlisHSPlasH/PBF/SimulationDataPBF.cpp b/SPlisHSPlasH/PBF/SimulationDataPBF.cpp
index c297103..5d6baaa 100644
--- a/SPlisHSPlasH/PBF/SimulationDataPBF.cpp
+++ b/SPlisHSPlasH/PBF/SimulationDataPBF.cpp
@@ -67,8 +67,8 @@ void SimulationDataPBF::reset()
{
m_deltaX[i][j].setZero();
m_lambda[i][j] = 0.0;
- getLastPosition(i, j) = fm->getPosition(i);
- getOldPosition(i, j) = fm->getPosition(i);
+ getLastPosition(i, j) = fm->getPosition(j);
+ getOldPosition(i, j) = fm->getPosition(j);
}
}
}
diff --git a/SPlisHSPlasH/PBF/TimeStepPBF.cpp b/SPlisHSPlasH/PBF/TimeStepPBF.cpp
index c0327e0..e6f2b4f 100644
--- a/SPlisHSPlasH/PBF/TimeStepPBF.cpp
+++ b/SPlisHSPlasH/PBF/TimeStepPBF.cpp
@@ -23,10 +23,27 @@ TimeStepPBF::TimeStepPBF() :
m_simulationData.init();
m_counter = 0;
m_velocityUpdateMethod = 0;
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "lambda", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getLambda(fluidModelIndex, i); } });
+ model->addField({ "deltaX", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDeltaX(fluidModelIndex, i)[0]; } });
+ }
}
TimeStepPBF::~TimeStepPBF(void)
{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("lambda");
+ model->removeFieldByName("deltaX");
+ }
}
void TimeStepPBF::initParameters()
@@ -153,7 +170,7 @@ void TimeStepPBF::pressureSolve()
for (unsigned int i = 0; i < nFluids; i++)
{
FluidModel *model = sim->getFluidModel(i);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
avg_density_err = 0.0;
pressureSolveIteration(i, avg_density_err);
@@ -178,7 +195,7 @@ void TimeStepPBF::pressureSolveIteration(const unsigned int fluidModelIndex, Rea
const unsigned int nFluids = sim->numberOfFluidModels();
const Real eps = 1.0e-6;
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
#pragma omp parallel default(shared)
{
diff --git a/SPlisHSPlasH/PCISPH/SimulationDataPCISPH.cpp b/SPlisHSPlasH/PCISPH/SimulationDataPCISPH.cpp
index c95f348..4e6377a 100644
--- a/SPlisHSPlasH/PCISPH/SimulationDataPCISPH.cpp
+++ b/SPlisHSPlasH/PCISPH/SimulationDataPCISPH.cpp
@@ -46,7 +46,7 @@ void SimulationDataPCISPH::init()
// Find prototype particle
// => particle with max. fluid neighbors
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
unsigned int index = 0;
unsigned int maxNeighbors = 0;
diff --git a/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp b/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
index 90f485f..961df08 100644
--- a/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
+++ b/SPlisHSPlasH/PCISPH/TimeStepPCISPH.cpp
@@ -14,10 +14,29 @@ TimeStepPCISPH::TimeStepPCISPH() :
{
m_simulationData.init();
m_counter = 0;
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "pressure", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressure(fluidModelIndex, i); } });
+ model->addField({ "advected density", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getDensityAdv(fluidModelIndex, i); } });
+ model->addField({ "pressure acceleration", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressureAccel(fluidModelIndex, i)[0]; } });
+ }
}
TimeStepPCISPH::~TimeStepPCISPH(void)
{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("pressure");
+ model->removeFieldByName("advected density");
+ model->removeFieldByName("pressure acceleration");
+ }
}
void TimeStepPCISPH::step()
@@ -101,7 +120,7 @@ void TimeStepPCISPH::pressureSolve()
for (unsigned int i = 0; i < nFluids; i++)
{
FluidModel *model = sim->getFluidModel(i);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
avg_density_err = 0.0;
pressureSolveIteration(i, avg_density_err);
@@ -120,7 +139,7 @@ void TimeStepPCISPH::pressureSolveIteration(const unsigned int fluidModelIndex,
Simulation *sim = Simulation::getCurrent();
FluidModel *model = sim->getFluidModel(fluidModelIndex);
const int numParticles = (int)model->numActiveParticles();
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
const Real h2 = h*h;
const Real invH2 = static_cast<Real>(1.0) / h2;
diff --git a/SPlisHSPlasH/Simulation.cpp b/SPlisHSPlasH/Simulation.cpp
index 2d2df9f..728e617 100644
--- a/SPlisHSPlasH/Simulation.cpp
+++ b/SPlisHSPlasH/Simulation.cpp
@@ -63,7 +63,7 @@ Simulation::Simulation ()
m_neighborhoodSearch = nullptr;
m_timeStep = nullptr;
- m_simulationMethod = SimulationMethods::WCSPH;
+ m_simulationMethod = SimulationMethods::NumSimulationMethods;
m_simulationMethodChanged = NULL;
m_sim2D = false;
@@ -113,7 +113,8 @@ void Simulation::init(const Real particleRadius, const bool sim2D)
// init kernel
setParticleRadius(particleRadius);
- setSimulationMethod(static_cast<int>(SimulationMethods::DFSPH));
+ setValue(Simulation::KERNEL_METHOD, Simulation::ENUM_KERNEL_PRECOMPUTED_CUBIC);
+ setValue(Simulation::GRAD_KERNEL_METHOD, Simulation::ENUM_GRADKERNEL_PRECOMPUTED_CUBIC);
// Initialize neighborhood search
if (m_neighborhoodSearch == NULL)
@@ -316,7 +317,7 @@ void Simulation::setKernel(int val)
m_kernelFct = WendlandQuinticC2Kernel2D::W;
}
}
- updateBoundaryPsi();
+ updateBoundaryVolume();
}
void Simulation::updateTimeStepSize()
@@ -396,6 +397,7 @@ void Simulation::computeNonPressureForces()
fm->computeViscosity();
fm->computeVorticity();
fm->computeDragForce();
+ fm->computeElasticity();
}
STOP_TIMING_AVG
}
@@ -409,7 +411,7 @@ void Simulation::reset()
// reset boundary models
for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
getBoundaryModel(i)->reset();
- updateBoundaryPsi();
+ updateBoundaryVolume();
if (m_timeStep)
m_timeStep->reset();
@@ -543,14 +545,11 @@ void Simulation::addFluidModel(const std::string &id, const unsigned int nFluidP
}
-void Simulation::updateBoundaryPsi()
+void Simulation::updateBoundaryVolume()
{
if (m_neighborhoodSearch == nullptr)
return;
- // ToDo
- const Real density0 = 1000.0;
-
Simulation *sim = Simulation::getCurrent();
const unsigned int nFluids = sim->numberOfFluidModels();
@@ -562,7 +561,7 @@ void Simulation::updateBoundaryPsi()
// Search boundary neighborhood
// Activate only static boundaries
- LOG_INFO << "Initialize boundary psi";
+ LOG_INFO << "Initialize boundary volume";
m_neighborhoodSearch->set_active(false);
for (unsigned int i = 0; i < numberOfBoundaryModels(); i++)
{
@@ -576,7 +575,7 @@ void Simulation::updateBoundaryPsi()
for (unsigned int body = 0; body < numberOfBoundaryModels(); body++)
{
if (!getBoundaryModel(body)->getRigidBodyObject()->isDynamic())
- getBoundaryModel(body)->computeBoundaryPsi(density0);
+ getBoundaryModel(body)->computeBoundaryVolume();
}
//////////////////////////////////////////////////////////////////////////
@@ -592,7 +591,7 @@ void Simulation::updateBoundaryPsi()
{
m_neighborhoodSearch->set_active(body + nFluids, true, true);
m_neighborhoodSearch->find_neighbors();
- getBoundaryModel(body)->computeBoundaryPsi(density0);
+ getBoundaryModel(body)->computeBoundaryVolume();
}
}
@@ -605,4 +604,4 @@ void Simulation::updateBoundaryPsi()
for (unsigned int j = numberOfFluidModels(); j < m_neighborhoodSearch->point_sets().size(); j++)
m_neighborhoodSearch->set_active(i, j, true);
}
-}
\ No newline at end of file
+}
diff --git a/SPlisHSPlasH/Simulation.h b/SPlisHSPlasH/Simulation.h
index f535b52..6045692 100644
--- a/SPlisHSPlasH/Simulation.h
+++ b/SPlisHSPlasH/Simulation.h
@@ -145,7 +145,7 @@ namespace SPH
BoundaryModel *getBoundaryModel(const unsigned int index) { return m_boundaryModels[index]; }
BoundaryModel *getBoundaryModelFromPointSet(const unsigned int pointSetIndex) { return static_cast<BoundaryModel*>(m_neighborhoodSearch->point_set(pointSetIndex).get_user_data()); }
const unsigned int numberOfBoundaryModels() const { return static_cast<unsigned int>(m_boundaryModels.size()); }
- void updateBoundaryPsi();
+ void updateBoundaryVolume();
int getKernel() const { return m_kernelMethod; }
void setKernel(int val);
diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp
index 78a8649..cc28fec 100644
--- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp
+++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Akinci2013.cpp
@@ -8,10 +8,13 @@ SurfaceTension_Akinci2013::SurfaceTension_Akinci2013(FluidModel *model) :
SurfaceTensionBase(model)
{
m_normals.resize(model->numParticles(), Vector3r::Zero());
+
+ model->addField({ "normal", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_normals[i][0]; } });
}
SurfaceTension_Akinci2013::~SurfaceTension_Akinci2013(void)
{
+ m_model->removeFieldByName("normal");
m_normals.clear();
}
@@ -51,7 +54,7 @@ void SurfaceTension_Akinci2013::computeNormals()
void SurfaceTension_Akinci2013::step()
{
Simulation *sim = Simulation::getCurrent();
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const Real supportRadius = sim->getSupportRadius();
const unsigned int numParticles = m_model->numActiveParticles();
const Real k = m_surfaceTension;
@@ -108,7 +111,7 @@ void SurfaceTension_Akinci2013::step()
if (length2 > 1.0e-9)
{
xixj = ((Real) 1.0 / sqrt(length2)) * xixj;
- ai -= k * bm_neighbor->getBoundaryPsi(neighborIndex) * xixj * AdhesionKernel::W(xi - xj);
+ ai -= k * density0 * bm_neighbor->getVolume(neighborIndex) * xixj * AdhesionKernel::W(xi - xj);
}
)
}
diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp
index 88e9b8f..40614bb 100644
--- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp
+++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_Becker2007.cpp
@@ -24,6 +24,7 @@ void SurfaceTension_Becker2007::step()
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
+ const Real density0 = model->getDensity0();
// Compute forces
#pragma omp parallel default(shared)
@@ -53,9 +54,9 @@ void SurfaceTension_Becker2007::step()
const Vector3r xixj = xi - xj;
const Real r2 = xixj.dot(xixj);
if (r2 > diameter2)
- ai -= k / m_model->getMass(i) * bm_neighbor->getBoundaryPsi(neighborIndex) * (xi - xj) * sim->W(xi - xj);
+ ai -= k / m_model->getMass(i) * density0 * bm_neighbor->getVolume(neighborIndex) * (xi - xj) * sim->W(xi - xj);
else
- ai -= k / m_model->getMass(i) * bm_neighbor->getBoundaryPsi(neighborIndex) * (xi - xj) * sim->W(Vector3r(diameter, 0.0, 0.0));
+ ai -= k / m_model->getMass(i) * density0 * bm_neighbor->getVolume(neighborIndex) * (xi - xj) * sim->W(Vector3r(diameter, 0.0, 0.0));
)
}
}
diff --git a/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp b/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp
index 628b337..cb36373 100644
--- a/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp
+++ b/SPlisHSPlasH/SurfaceTension/SurfaceTension_He2014.cpp
@@ -9,10 +9,16 @@ SurfaceTension_He2014::SurfaceTension_He2014(FluidModel *model) :
{
m_color.resize(model->numParticles(), 0.0);
m_gradC2.resize(model->numParticles(), 0.0);
+
+ model->addField({ "color", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &m_color[i]; } });
+ model->addField({ "gradC2", FieldType::Scalar, [&](const unsigned int i) -> Real* { return &m_gradC2[i]; } });
}
SurfaceTension_He2014::~SurfaceTension_He2014(void)
{
+ m_model->removeFieldByName("color");
+ m_model->removeFieldByName("gradC2");
+
m_color.clear();
m_gradC2.clear();
}
@@ -23,7 +29,7 @@ void SurfaceTension_He2014::step()
Simulation *sim = Simulation::getCurrent();
const unsigned int numParticles = m_model->numActiveParticles();
const Real k = m_surfaceTension;
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
@@ -50,7 +56,7 @@ void SurfaceTension_He2014::step()
// Boundary
//////////////////////////////////////////////////////////////////////////
forall_boundary_neighbors(
- ci += bm_neighbor->getBoundaryPsi(neighborIndex) / density0 * sim->W(xi - xj);
+ ci += bm_neighbor->getVolume(neighborIndex) * sim->W(xi - xj);
)
}
}
@@ -104,7 +110,7 @@ void SurfaceTension_He2014::step()
// Boundary
//////////////////////////////////////////////////////////////////////////
forall_boundary_neighbors(
- ai += factor*bm_neighbor->getBoundaryPsi(neighborIndex) / density0 * gradC2_i * sim->gradW(xi - xj);
+ ai += factor* bm_neighbor->getVolume(neighborIndex) * gradC2_i * sim->gradW(xi - xj);
)
}
}
diff --git a/SPlisHSPlasH/Utilities/MathFunctions.cpp b/SPlisHSPlasH/Utilities/MathFunctions.cpp
new file mode 100644
index 0000000..273d770
--- /dev/null
+++ b/SPlisHSPlasH/Utilities/MathFunctions.cpp
@@ -0,0 +1,231 @@
+#include "MathFunctions.h"
+#include <cfloat>
+
+using namespace SPH;
+
+
+// ----------------------------------------------------------------------------------------------
+void MathFunctions::extractRotation(const Matrix3r &A, Quaternionr &q, const unsigned int maxIter)
+{
+ for (unsigned int iter = 0; iter < maxIter; iter++)
+ {
+ Matrix3r R = q.matrix();
+ Vector3r omega = (R.col(0).cross(A.col(0)) + R.col(1).cross(A.col(1)) + R.col(2).cross(A.col(2))) *
+ (1.0 / fabs(R.col(0).dot(A.col(0)) + R.col(1).dot(A.col(1)) + R.col(2).dot(A.col(2)) + 1.0e-9));
+ Real w = omega.norm();
+ if (w < 1.0e-9)
+ break;
+ q = Quaternionr(AngleAxisr(w, (1.0 / w) * omega)) * q;
+ q.normalize();
+ }
+}
+
+void MathFunctions::pseudoInverse(const Matrix3r &a, Matrix3r &res)
+{
+ const Real epsilon = std::numeric_limits<Real>::epsilon();
+ const Eigen::JacobiSVD<Matrix3r> svd(a, Eigen::ComputeFullU | Eigen::ComputeFullV);
+ const Real tolerance = epsilon * std::max(a.cols(), a.rows()) * svd.singularValues().array().abs()(0);
+ res = svd.matrixV() * (svd.singularValues().array().abs() > tolerance).select(svd.singularValues().array().inverse(), 0).matrix().asDiagonal() * svd.matrixU().adjoint();
+}
+
+
+/** Perform a singular value decomposition of matrix A: A = U * sigma * V^T.
+* This function returns two proper rotation matrices U and V^T which do not
+* contain a reflection. Reflections are corrected by the inversion handling
+* proposed by Irving et al. 2004.
+*/
+void MathFunctions::svdWithInversionHandling(const Matrix3r &A, Vector3r &sigma, Matrix3r &U, Matrix3r &VT)
+{
+
+ Matrix3r AT_A, V;
+ AT_A = A.transpose() * A;
+
+ Vector3r S;
+
+ // Eigen decomposition of A^T * A
+ eigenDecomposition(AT_A, V, S);
+
+ // Detect if V is a reflection .
+ // Make a rotation out of it by multiplying one column with -1.
+ const Real detV = V.determinant();
+ if (detV < 0.0)
+ {
+ Real minLambda = REAL_MAX;
+ unsigned char pos = 0;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ if (S[l] < minLambda)
+ {
+ pos = l;
+ minLambda = S[l];
+ }
+ }
+ V(0, pos) = -V(0, pos);
+ V(1, pos) = -V(1, pos);
+ V(2, pos) = -V(2, pos);
+ }
+
+ if (S[0] < 0.0) S[0] = 0.0; // safety for sqrt
+ if (S[1] < 0.0) S[1] = 0.0;
+ if (S[2] < 0.0) S[2] = 0.0;
+
+ sigma[0] = sqrt(S[0]);
+ sigma[1] = sqrt(S[1]);
+ sigma[2] = sqrt(S[2]);
+
+ VT = V.transpose();
+
+ //
+ // Check for values of hatF near zero
+ //
+ unsigned char chk = 0;
+ unsigned char pos = 0;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ if (fabs(sigma[l]) < 1.0e-4)
+ {
+ pos = l;
+ chk++;
+ }
+ }
+
+ if (chk > 0)
+ {
+ if (chk > 1)
+ {
+ U.setIdentity();
+ }
+ else
+ {
+ U = A * V;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ if (l != pos)
+ {
+ for (unsigned char m = 0; m < 3; m++)
+ {
+ U(m, l) *= static_cast<Real>(1.0) / sigma[l];
+ }
+ }
+ }
+
+ Vector3r v[2];
+ unsigned char index = 0;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ if (l != pos)
+ {
+ v[index++] = Vector3r(U(0, l), U(1, l), U(2, l));
+ }
+ }
+ Vector3r vec = v[0].cross(v[1]);
+ vec.normalize();
+ U(0, pos) = vec[0];
+ U(1, pos) = vec[1];
+ U(2, pos) = vec[2];
+ }
+ }
+ else
+ {
+ Vector3r sigmaInv(static_cast<Real>(1.0) / sigma[0], static_cast<Real>(1.0) / sigma[1], static_cast<Real>(1.0) / sigma[2]);
+ U = A * V;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ for (unsigned char m = 0; m < 3; m++)
+ {
+ U(m, l) *= sigmaInv[l];
+ }
+ }
+ }
+
+ const Real detU = U.determinant();
+
+ // U is a reflection => inversion
+ if (detU < 0.0)
+ {
+ //std::cout << "Inversion!\n";
+ Real minLambda = REAL_MAX;
+ unsigned char pos = 0;
+ for (unsigned char l = 0; l < 3; l++)
+ {
+ if (sigma[l] < minLambda)
+ {
+ pos = l;
+ minLambda = sigma[l];
+ }
+ }
+
+ // invert values of smallest singular value
+ sigma[pos] = -sigma[pos];
+ U(0, pos) = -U(0, pos);
+ U(1, pos) = -U(1, pos);
+ U(2, pos) = -U(2, pos);
+ }
+}
+
+// ----------------------------------------------------------------------------------------------
+void MathFunctions::eigenDecomposition(const Matrix3r &A, Matrix3r &eigenVecs, Vector3r &eigenVals)
+{
+ const int numJacobiIterations = 10;
+ const Real epsilon = static_cast<Real>(1e-15);
+
+ Matrix3r D = A;
+
+ // only for symmetric matrices!
+ eigenVecs.setIdentity(); // unit matrix
+ int iter = 0;
+ while (iter < numJacobiIterations) { // 3 off diagonal elements
+ // find off diagonal element with maximum modulus
+ int p, q;
+ Real a, max;
+ max = fabs(D(0, 1));
+ p = 0; q = 1;
+ a = fabs(D(0, 2));
+ if (a > max) { p = 0; q = 2; max = a; }
+ a = fabs(D(1, 2));
+ if (a > max) { p = 1; q = 2; max = a; }
+ // all small enough -> done
+ if (max < epsilon) break;
+ // rotate matrix with respect to that element
+ jacobiRotate(D, eigenVecs, p, q);
+ iter++;
+ }
+ eigenVals[0] = D(0, 0);
+ eigenVals[1] = D(1, 1);
+ eigenVals[2] = D(2, 2);
+}
+
+// ----------------------------------------------------------------------------------------------
+void MathFunctions::jacobiRotate(Matrix3r &A, Matrix3r &R, int p, int q)
+{
+ // rotates A through phi in pq-plane to set A(p,q) = 0
+ // rotation stored in R whose columns are eigenvectors of A
+ if (A(p, q) == 0.0)
+ return;
+
+ Real d = (A(p, p) - A(q, q)) / (static_cast<Real>(2.0)*A(p, q));
+ Real t = static_cast<Real>(1.0) / (fabs(d) + sqrt(d*d + static_cast<Real>(1.0)));
+ if (d < 0.0) t = -t;
+ Real c = static_cast<Real>(1.0) / sqrt(t*t + 1);
+ Real s = t*c;
+ A(p, p) += t*A(p, q);
+ A(q, q) -= t*A(p, q);
+ A(p, q) = A(q, p) = 0.0;
+ // transform A
+ int k;
+ for (k = 0; k < 3; k++) {
+ if (k != p && k != q) {
+ Real Akp = c*A(k, p) + s*A(k, q);
+ Real Akq = -s*A(k, p) + c*A(k, q);
+ A(k, p) = A(p, k) = Akp;
+ A(k, q) = A(q, k) = Akq;
+ }
+ }
+ // store rotation in R
+ for (k = 0; k < 3; k++) {
+ Real Rkp = c*R(k, p) + s*R(k, q);
+ Real Rkq = -s*R(k, p) + c*R(k, q);
+ R(k, p) = Rkp;
+ R(k, q) = Rkq;
+ }
+}
diff --git a/SPlisHSPlasH/Utilities/MathFunctions.h b/SPlisHSPlasH/Utilities/MathFunctions.h
new file mode 100644
index 0000000..1491d0e
--- /dev/null
+++ b/SPlisHSPlasH/Utilities/MathFunctions.h
@@ -0,0 +1,26 @@
+#ifndef MATH_FUNCTIONS_H
+#define MATH_FUNCTIONS_H
+
+#include "SPlisHSPlasH/Common.h"
+
+// ------------------------------------------------------------------------------------
+namespace SPH
+{
+ class MathFunctions
+ {
+ public:
+ /** Implementation of the paper: \n
+ * Matthias M�ller, Jan Bender, Nuttapong Chentanez and Miles Macklin,
+ * "A Robust Method to Extract the Rotational Part of Deformations",
+ * ACM SIGGRAPH Motion in Games, 2016
+ */
+ static void extractRotation(const Matrix3r &A, Quaternionr &q, const unsigned int maxIter);
+
+ static void pseudoInverse(const Matrix3r &a, Matrix3r &res);
+ static void svdWithInversionHandling(const Matrix3r &A, Vector3r &sigma, Matrix3r &U, Matrix3r &VT);
+ static void eigenDecomposition(const Matrix3r &A, Matrix3r &eigenVecs, Vector3r &eigenVals);
+ static void jacobiRotate(Matrix3r &A, Matrix3r &R, int p, int q);
+ };
+}
+
+#endif
\ No newline at end of file
diff --git a/SPlisHSPlasH/Utilities/SceneLoader.cpp b/SPlisHSPlasH/Utilities/SceneLoader.cpp
index f804ebd..5749952 100644
--- a/SPlisHSPlasH/Utilities/SceneLoader.cpp
+++ b/SPlisHSPlasH/Utilities/SceneLoader.cpp
@@ -328,6 +328,35 @@ void SceneLoader::readParameterObject(const std::string &key, ParameterObject *p
static_cast<VectorParameter<Real>*>(paramBase)->setValue(val.data());
}
}
+ else if (paramBase->getType() == ParameterBase::STRING)
+ {
+ std::string val;
+ if (readValue(config[paramBase->getName()], val))
+ static_cast<StringParameter*>(paramBase)->setValue(val);
+ }
}
}
}
+
+Utilities::SceneLoader::ColoringData Utilities::SceneLoader::readColoringInfo(const std::string &key)
+{
+ ColoringData data;
+ data.colorField = "velocity";
+ data.colorMapType = 0;
+ data.minVal = 0.0;
+ data.maxVal = 5.0;
+
+ //////////////////////////////////////////////////////////////////////////
+ // read configuration
+ //////////////////////////////////////////////////////////////////////////
+ if (m_jsonData.find(key) != m_jsonData.end())
+ {
+ nlohmann::json config = m_jsonData[key];
+#
+ readValue(config["renderMinValue"], data.minVal);
+ readValue(config["renderMaxValue"], data.maxVal);
+ readValue(config["colorField"], data.colorField);
+ readValue(config["colorMapType"], data.colorMapType);
+ }
+ return data;
+}
diff --git a/SPlisHSPlasH/Utilities/SceneLoader.h b/SPlisHSPlasH/Utilities/SceneLoader.h
index 6b33483..ee46b84 100644
--- a/SPlisHSPlasH/Utilities/SceneLoader.h
+++ b/SPlisHSPlasH/Utilities/SceneLoader.h
@@ -82,6 +82,15 @@ namespace Utilities
Real timeStepSize;
};
+ /** \brief Struct to store particle coloring information */
+ struct ColoringData
+ {
+ std::string colorField;
+ unsigned int colorMapType;
+ Real minVal;
+ Real maxVal;
+ };
+
void readScene(const char *fileName, Scene &scene);
@@ -148,6 +157,7 @@ namespace Utilities
}
void readParameterObject(const std::string &key, GenParam::ParameterObject *paramObj);
+ ColoringData readColoringInfo(const std::string &key);
};
template <>
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp
index 7a1287c..d6934ca 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Bender2017.cpp
@@ -21,10 +21,18 @@ Viscosity_Bender2017::Viscosity_Bender2017(FluidModel *model) :
m_iterations = 0;
m_maxIter = 50;
m_maxError = 0.01;
+
+ model->addField({ "target strain rate", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_targetStrainRate[i][0]; } });
+ model->addField({ "viscosity factor", FieldType::Matrix6, [&](const unsigned int i) -> Real* { return &m_viscosityFactor[i](0,0); } });
+ model->addField({ "viscosity lambda", FieldType::Vector6, [&](const unsigned int i) -> Real* { return &m_viscosityLambda[i][0]; } });
}
Viscosity_Bender2017::~Viscosity_Bender2017(void)
{
+ m_model->removeFieldByName("target strain rate");
+ m_model->removeFieldByName("viscosity factor");
+ m_model->removeFieldByName("viscosity lambda");
+
m_targetStrainRate.clear();
m_viscosityFactor.clear();
m_viscosityLambda.clear();
@@ -61,6 +69,7 @@ void Viscosity_Bender2017::step()
const unsigned int maxIter = m_maxIter;
const Real maxError = m_maxError;
const Real maxError2 = maxError*maxError;
+ const Real density0 = model->getDensity0();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
@@ -187,7 +196,7 @@ void Viscosity_Bender2017::step()
//////////////////////////////////////////////////////////////////////////
forall_boundary_neighbors(
const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
- ai -= invH * 0.1 * m_viscosity * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) * (vi - vj)* sim->W(xi - xj);
+ ai -= invH * 0.1 * m_viscosity * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi - vj)* sim->W(xi - xj);
)
}
}
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp
index 0409c6b..9e1a8fc 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Peer2015.cpp
@@ -20,10 +20,14 @@ Viscosity_Peer2015::Viscosity_Peer2015(FluidModel *model) :
m_iterations = 0;
m_maxIter = 50;
m_maxError = 0.01;
+
+ model->addField({ "target nablaV", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_targetNablaV[i](0,0); } });
}
Viscosity_Peer2015::~Viscosity_Peer2015(void)
{
+ m_model->removeFieldByName("target nablaV");
+
m_targetNablaV.clear();
}
@@ -88,7 +92,7 @@ void Viscosity_Peer2015::step()
Simulation *sim = Simulation::getCurrent();
const int numParticles = (int) m_model->numActiveParticles();
const Real viscosity = static_cast<Real>(1.0) - m_viscosity;
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp
index 98adad2..9971302 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Peer2016.cpp
@@ -27,10 +27,16 @@ Viscosity_Peer2016::Viscosity_Peer2016(FluidModel *model) :
m_maxErrorV = 0.01;
m_maxIterOmega = 50;
m_maxErrorOmega = 0.01;
+
+ model->addField({ "target nablaV", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_targetNablaV[i](0,0); } });
+ model->addField({ "omega (visco)", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_omega[i][0]; } });
}
Viscosity_Peer2016::~Viscosity_Peer2016(void)
{
+ m_model->removeFieldByName("target nablaV");
+ m_model->removeFieldByName("omega (visco)");
+
m_targetNablaV.clear();
m_omega.clear();
}
@@ -174,7 +180,7 @@ void Viscosity_Peer2016::step()
Simulation *sim = Simulation::getCurrent();
const int numParticles = (int) m_model->numActiveParticles();
const Real viscosity = static_cast<Real>(1.0) - m_viscosity;
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp
index d37b411..7297e4f 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Standard.cpp
@@ -21,6 +21,9 @@ void Viscosity_Standard::step()
const Real h2 = h*h;
const unsigned int nFluids = sim->numberOfFluidModels();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ Real d = 10.0;
+ if (sim->is2DSimulation())
+ d = 6.0;
#pragma omp parallel default(shared)
{
@@ -41,7 +44,7 @@ void Viscosity_Standard::step()
// Viscosity
const Real density_j = fm_neighbor->getDensity(neighborIndex);
const Vector3r xixj = xi - xj;
- ai += 10.0 * m_viscosity * (fm_neighbor->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * sim->gradW(xi - xj);
+ ai += d * m_viscosity * (fm_neighbor->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * sim->gradW(xi - xj);
)
////////////////////////////////////////////////////////////////////////////
@@ -56,7 +59,7 @@ void Viscosity_Standard::step()
// const Vector3r &xj = bm_neighbor->getPosition(neighborIndex);
// const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
// const Vector3r xixj = xi - xj;
- // ai += 10.0 * viscosity * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) * (vi) * (xixj.dot(sim->gradW(xi - xj))) / (xixj.squaredNorm() + 0.01*h2);
+ // ai += d * viscosity * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi) * (xixj.dot(sim->gradW(xi - xj))) / (xixj.squaredNorm() + 0.01*h2);
// }
//}
}
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp
index 2c51518..a682464 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Takahashi2015.cpp
@@ -20,10 +20,16 @@ Viscosity_Takahashi2015::Viscosity_Takahashi2015(FluidModel *model) :
m_maxIter = 100;
m_maxError = 0.01;
m_iterations = 0;
+
+ model->addField({ "viscous stress", FieldType::Matrix3, [&](const unsigned int i) -> Real* { return &m_viscousStress[i](0,0); } });
+ model->addField({ "accel (visco)", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_accel[i][0]; } });
}
Viscosity_Takahashi2015::~Viscosity_Takahashi2015(void)
{
+ m_model->removeFieldByName("viscous stress");
+ m_model->removeFieldByName("accel (visco)");
+
m_viscousStress.clear();
m_accel.clear();
}
@@ -132,7 +138,7 @@ void Viscosity_Takahashi2015::step()
{
Simulation *sim = Simulation::getCurrent();
const int numParticles = (int) m_model->numActiveParticles();
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
const unsigned int nFluids = sim->numberOfFluidModels();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
@@ -208,7 +214,7 @@ void Viscosity_Takahashi2015::step()
const unsigned int neighborIndex = sim->getNeighbor(fluidModelIndex, pid, i, j);
const Vector3r &xj = bm_neighbor->getPosition(neighborIndex);
const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
- ai -= invH * 0.1 * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) * (vi - vj)* sim->W(xi - xj);
+ ai -= invH * 0.1 * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi - vj)* sim->W(xi - xj);
}
}
}
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp b/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp
index c20d28b..620a1ad 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_Weiler2018.cpp
@@ -21,10 +21,14 @@ Viscosity_Weiler2018::Viscosity_Weiler2018(FluidModel *model) :
m_boundaryViscosity = 0.0;
m_vDiff.resize(model->numParticles(), Vector3r::Zero());
+
+ model->addField({ "velocity difference", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_vDiff[i][0]; } });
}
Viscosity_Weiler2018::~Viscosity_Weiler2018(void)
{
+ m_model->removeFieldByName("velocity difference");
+
m_vDiff.clear();
}
@@ -69,6 +73,11 @@ void Viscosity_Weiler2018::matrixVecProd(const Real* vec, Real *result, void *us
const Real dt = TimeManager::getCurrent()->getTimeStepSize();
const Real mu = visco->m_viscosity;
const Real mub = visco->m_boundaryViscosity;
+ const Real density0 = model->getDensity0();
+
+ Real d = 10.0;
+ if (sim->is2DSimulation())
+ d = 6.0;
#pragma omp parallel default(shared)
{
@@ -91,7 +100,7 @@ void Viscosity_Weiler2018::matrixVecProd(const Real* vec, Real *result, void *us
const Vector3r &vj = Eigen::Map<const Vector3r>(&vec[3 * neighborIndex]);
const Vector3r xixj = xi - xj;
- ai += 10.0 * mu * (model->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW;
+ ai += d * mu * (model->getMass(neighborIndex) / density_j) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW;
)
//////////////////////////////////////////////////////////////////////////
@@ -102,7 +111,7 @@ void Viscosity_Weiler2018::matrixVecProd(const Real* vec, Real *result, void *us
const Vector3r gradW = sim->gradW(xi - xj);
const Vector3r xixj = xi - xj;
- ai += 10.0 * mub * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW;
+ ai += d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * (vi - vj).dot(xixj) / (xixj.squaredNorm() + 0.01*h2) * gradW;
)
result[3 * i] = vec[3 * i] - dt / density_i*ai[0];
@@ -121,7 +130,11 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Matrix3r
FluidModel *model = visco->getModel();
const unsigned int nFluids = sim->numberOfFluidModels();
const unsigned int fluidModelIndex = model->getPointSetIndex();
+ const Real density0 = model->getDensity0();
+ Real d = 10.0;
+ if (sim->is2DSimulation())
+ d = 6.0;
const Real h = sim->getSupportRadius();
const Real h2 = h*h;
@@ -142,7 +155,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Matrix3r
const Real density_j = model->getDensity(neighborIndex);
const Vector3r gradW = sim->gradW(xi - xj);
const Vector3r xixj = xi - xj;
- result += 10.0 * mu * (model->getMass(neighborIndex) / density_j) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
+ result += d * mu * (model->getMass(neighborIndex) / density_j) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
)
//////////////////////////////////////////////////////////////////////////
@@ -153,7 +166,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Matrix3r
const Vector3r gradW = sim->gradW(xi - xj);
const Vector3r xixj = xi - xj;
- result += 10.0 * mub * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
+ result += d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
)
result = Matrix3r::Identity() - (dt / density_i) * result;
}
@@ -173,7 +186,11 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Vector3r
const Real mub = visco->m_boundaryViscosity;
const unsigned int nFluids = sim->numberOfFluidModels();
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
+ const Real density0 = model->getDensity0();
+ Real d = 10.0;
+ if (sim->is2DSimulation())
+ d = 6.0;
const Real density_i = model->getDensity(i);
@@ -188,7 +205,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Vector3r
const Real density_j = model->getDensity(neighborIndex);
const Vector3r gradW = model->gradW(xi - xj);
const Vector3r xixj = xi - xj;
- Matrix3r r = 10.0 * mu * (model->getMass(neighborIndex) / density_j) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
+ Matrix3r r = d * mu * (model->getMass(neighborIndex) / density_j) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
result += r.diagonal();
)
@@ -200,7 +217,7 @@ void Viscosity_Weiler2018::diagonalMatrixElement(const unsigned int i, Vector3r
const Vector3r gradW = model->gradW(xi - xj);
const Vector3r xixj = xi - xj;
- Matrix3r r = 10.0 * mub * (model->getBoundaryPsi(pid, neighborIndex) / density_i) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
+ Matrix3r r = d * mub * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) / (xixj.squaredNorm() + 0.01*h2) * (gradW * xixj.transpose());
result += r.diagonal();
)
result = Vector3r::Ones() - (dt / density_i) * result;
@@ -215,7 +232,7 @@ void Viscosity_Weiler2018::step()
// prevent solver from running with a zero-length vector
if (numParticles == 0)
return;
- const Real density0 = m_model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = m_model->getDensity0();
const Real h = TimeManager::getCurrent()->getTimeStepSize();
//////////////////////////////////////////////////////////////////////////
diff --git a/SPlisHSPlasH/Viscosity/Viscosity_XSPH.cpp b/SPlisHSPlasH/Viscosity/Viscosity_XSPH.cpp
index 4fd861d..890dec0 100644
--- a/SPlisHSPlasH/Viscosity/Viscosity_XSPH.cpp
+++ b/SPlisHSPlasH/Viscosity/Viscosity_XSPH.cpp
@@ -63,7 +63,7 @@ void Viscosity_XSPH::step()
// const unsigned int neighborIndex = sim->getNeighbor(pid, i, j);
// const Vector3r &xj = bm_neighbor->getPosition(neighborIndex);
// const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
- // ai -= invH * viscosity * (bm_neighbor->getBoundaryPsi(neighborIndex) / density_i) * (vi)* sim->W(xi - xj);
+ // ai -= invH * viscosity * bm_neighbor->getVolume(neighborIndex) * (vi)* sim->W(xi - xj);
// }
//}
}
diff --git a/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp b/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp
index cabfa16..2ed123c 100644
--- a/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp
+++ b/SPlisHSPlasH/Vorticity/MicropolarModel_Bender2017.cpp
@@ -17,10 +17,14 @@ MicropolarModel_Bender2017::MicropolarModel_Bender2017(FluidModel *model) :
m_angularAcceleration.resize(model->numParticles(), Vector3r::Zero());
m_inertiaInverse = 0.5;
m_viscosityOmega = 0.1;
+
+ model->addField({ "angular velocity", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_omega[i][0]; } });
}
MicropolarModel_Bender2017::~MicropolarModel_Bender2017(void)
{
+ m_model->removeFieldByName("angular velocity");
+
m_omega.clear();
m_angularAcceleration.clear();
}
@@ -51,13 +55,21 @@ void MicropolarModel_Bender2017::step()
const unsigned int fluidModelIndex = m_model->getPointSetIndex();
const unsigned int nFluids = sim->numberOfFluidModels();
FluidModel *model = m_model;
+ const Real density0 = model->getDensity0();
- const Real h = TimeManager::getCurrent()->getTimeStepSize();
- const Real invH = static_cast<Real>(1.0) / h;
+ const Real dt = TimeManager::getCurrent()->getTimeStepSize();
+ const Real invDt = static_cast<Real>(1.0) / dt;
const Real nu_t = m_vorticityCoeff;
const Real zeta = m_viscosityOmega;
+ const Real h = sim->getSupportRadius();
+ const Real h2 = h*h;
+
+ Real d = 10.0;
+ if (sim->is2DSimulation())
+ d = 6.0;
+
#pragma omp parallel default(shared)
{
#pragma omp for schedule(static)
@@ -86,13 +98,47 @@ void MicropolarModel_Bender2017::step()
const Vector3r omegaij = omegai - omegaj;
const Vector3r gradW = sim->gradW(xij);
- // XSPH for angular velocity field
- angAcceli -= invH * m_inertiaInverse * zeta * (m_model->getMass(neighborIndex) / density_j) * omegaij * sim->W(xij);
-
- // symmetric curl
- ai -= nu_t * density_i * m_model->getMass(neighborIndex) * ((omegai / density_i2 + omegaj / density_j2).cross(gradW));
- angAcceli -= nu_t * density_i * m_inertiaInverse * (m_model->getMass(neighborIndex) * (vi / density_i2 + vj / density_j2).cross(gradW));
- )
+ // XSPH for angular velocity field
+ angAcceli -= invDt * m_inertiaInverse * zeta * (m_model->getMass(neighborIndex) / density_j) * omegaij * sim->W(xij);
+
+ //// Viscosity
+ //angAcceli += d * m_inertiaInverse * zeta * (m_model->getMass(neighborIndex) / density_i) * omegaij.dot(xij) / (xij.squaredNorm() + 0.01*h2) * gradW;
+
+ // difference curl
+ ai += nu_t * 1.0/density_i * m_model->getMass(neighborIndex) * (omegaij.cross(gradW));
+ angAcceli += nu_t * 1.0/density_i * m_inertiaInverse * (m_model->getMass(neighborIndex) * (vi - vj).cross(gradW));
+
+// // symmetric curl
+// ai -= nu_t * density_i * m_model->getMass(neighborIndex) * ((omegai / density_i2 + omegaj / density_j2).cross(gradW));
+// angAcceli -= nu_t * density_i * m_inertiaInverse * (m_model->getMass(neighborIndex) * (vi / density_i2 + vj / density_j2).cross(gradW));
+ );
+
+ //////////////////////////////////////////////////////////////////////////
+ // Boundary
+ //////////////////////////////////////////////////////////////////////////
+ forall_boundary_neighbors(
+ const Vector3r &vj = bm_neighbor->getVelocity(neighborIndex);
+ const Vector3r &omegaj = Vector3r::Zero();//m_omega[neighborIndex];
+
+ // Viscosity
+ const Vector3r xij = xi - xj;
+ const Vector3r omegaij = omegai - omegaj;
+ const Vector3r gradW = sim->gradW(xij);
+
+ // XSPH for angular velocity field
+ //angAcceli -= invDt * m_inertiaInverse * zeta * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * omegaij * sim->W(xij);
+
+ // Viscosity
+ //angAcceli += d * m_inertiaInverse * zeta * (density0 * bm_neighbor->getVolume(neighborIndex) / density_i) * omegaij.dot(xij) / (xij.squaredNorm() + 0.01*h2) * gradW;
+
+ // difference curl
+ ai += nu_t * 1.0 / density_i * density0 * bm_neighbor->getVolume(neighborIndex) * (omegaij.cross(gradW));
+ angAcceli += nu_t * 1.0 / density_i * m_inertiaInverse * (density0 * bm_neighbor->getVolume(neighborIndex) * (vi - vj).cross(gradW));
+
+// // symmetric curl
+// ai -= nu_t * density_i * density0 * bm_neighbor->getVolume(neighborIndex) * ((omegai / density_i2 + omegaj / density_i2).cross(gradW));
+// angAcceli -= nu_t * density_i * m_inertiaInverse * (density0 * bm_neighbor->getVolume(neighborIndex) * (vi / density_i2 + vj / density_i2).cross(gradW));
+ );
angAcceli -= 2.0 * m_inertiaInverse * nu_t * omegai;
}
}
@@ -102,7 +148,7 @@ void MicropolarModel_Bender2017::step()
#pragma omp for schedule(static)
for (int i = 0; i < (int)numParticles; i++)
{
- m_omega[i] += h*m_angularAcceleration[i];
+ m_omega[i] += dt*m_angularAcceleration[i];
}
}
}
diff --git a/SPlisHSPlasH/Vorticity/VorticityConfinement.cpp b/SPlisHSPlasH/Vorticity/VorticityConfinement.cpp
index 79df01d..5fcaccc 100644
--- a/SPlisHSPlasH/Vorticity/VorticityConfinement.cpp
+++ b/SPlisHSPlasH/Vorticity/VorticityConfinement.cpp
@@ -10,10 +10,14 @@ VorticityConfinement::VorticityConfinement(FluidModel *model) :
{
m_omega.resize(model->numParticles(), Vector3r::Zero());
m_normOmega.resize(model->numParticles(), 0.0);
+
+ model->addField({ "angular velocity", FieldType::Vector3, [&](const unsigned int i) -> Real* { return &m_omega[i][0]; } });
}
VorticityConfinement::~VorticityConfinement(void)
{
+ m_model->removeFieldByName("angular velocity");
+
m_omega.clear();
m_normOmega.clear();
}
diff --git a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
index 85e31b5..ebc16cd 100644
--- a/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
+++ b/SPlisHSPlasH/WCSPH/TimeStepWCSPH.cpp
@@ -22,10 +22,27 @@ TimeStepWCSPH::TimeStepWCSPH() :
m_stiffness = 50.0;
m_exponent = 7.0;
+
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->addField({ "pressure", FieldType::Scalar, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressure(fluidModelIndex, i); } });
+ model->addField({ "pressure acceleration", FieldType::Vector3, [this, fluidModelIndex](const unsigned int i) -> Real* { return &m_simulationData.getPressureAccel(fluidModelIndex, i)[0]; } });
+ }
}
TimeStepWCSPH::~TimeStepWCSPH(void)
{
+ Simulation *sim = Simulation::getCurrent();
+ const unsigned int nModels = sim->numberOfFluidModels();
+ for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
+ {
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
+ model->removeFieldByName("pressure");
+ model->removeFieldByName("pressure acceleration");
+ }
}
void TimeStepWCSPH::initParameters()
@@ -65,7 +82,7 @@ void TimeStepWCSPH::step()
for (unsigned int fluidModelIndex = 0; fluidModelIndex < nModels; fluidModelIndex++)
{
FluidModel *model = sim->getFluidModel(fluidModelIndex);
- const Real density0 = model->getValue<Real>(FluidModel::DENSITY0);
+ const Real density0 = model->getDensity0();
#pragma omp parallel default(shared)
{
#pragma omp for schedule(static)
diff --git a/Demos/CMakeLists.txt b/Simulators/CMakeLists.txt
similarity index 63%
rename from Demos/CMakeLists.txt
rename to Simulators/CMakeLists.txt
index e66ca18..6e30108 100644
--- a/Demos/CMakeLists.txt
+++ b/Simulators/CMakeLists.txt
@@ -1,5 +1,5 @@
include(${PROJECT_PATH}/Visualization/CMakeLists.txt)
add_definitions(-DPBD_DATA_PATH="../data")
-subdirs(DynamicBoundaryDemo StaticBoundaryDemo)
+subdirs(DynamicBoundarySimulator StaticBoundarySimulator)
diff --git a/Demos/Common/DemoBase.cpp b/Simulators/Common/SimulatorBase.cpp
similarity index 70%
rename from Demos/Common/DemoBase.cpp
rename to Simulators/Common/SimulatorBase.cpp
index 931b256..d4638eb 100644
--- a/Demos/Common/DemoBase.cpp
+++ b/Simulators/Common/SimulatorBase.cpp
@@ -1,4 +1,4 @@
-#include "DemoBase.h"
+#include "SimulatorBase.h"
#include "Visualization/MiniGL.h"
#include "SPlisHSPlasH/Utilities/SceneLoader.h"
#include "Utilities/FileSystem.h"
@@ -18,6 +18,7 @@
#include "Utilities/SystemInfo.h"
#include "Visualization/colormaps/colormap_jet.h"
#include "Visualization/colormaps/colormap_plasma.h"
+#include "extern/partio/src/lib/Partio.h"
INIT_LOGGING
INIT_TIMING
@@ -28,40 +29,30 @@ using namespace std;
using namespace GenParam;
using namespace Utilities;
-int DemoBase::PAUSE = -1;
-int DemoBase::PAUSE_AT = -1;
-int DemoBase::STOP_AT = -1;
-int DemoBase::NUM_STEPS_PER_RENDER = -1;
-int DemoBase::PARTIO_EXPORT = -1;
-int DemoBase::PARTIO_EXPORT_FPS = -1;
-int DemoBase::RENDER_MIN_VALUE = -1;
-int DemoBase::RENDER_MAX_VALUE = -1;
-int DemoBase::RENDER_WALLS = -1;
-int DemoBase::RENDER_COLOR_FIELD = -1;
-int DemoBase::RENDER_COLOR_MAP_TYPE = -1;
-int DemoBase::ENUM_WALLS_NONE = -1;
-int DemoBase::ENUM_WALLS_PARTICLES_ALL = -1;
-int DemoBase::ENUM_WALLS_PARTICLES_NO_WALLS = -1;
-int DemoBase::ENUM_WALLS_GEOMETRY_ALL = -1;
-int DemoBase::ENUM_WALLS_GEOMETRY_NO_WALLS = -1;
-int DemoBase::ENUM_RENDER_NONE = -1;
-int DemoBase::ENUM_RENDER_VELOCITY = -1;
-int DemoBase::ENUM_RENDER_ANGULAR_VELOCITY = -1;
-int DemoBase::ENUM_RENDER_DENSITY = -1;
-int DemoBase::ENUM_COLORMAP_NONE = -1;
-int DemoBase::ENUM_COLORMAP_JET = -1;
-int DemoBase::ENUM_COLORMAP_PLASMA = -1;
+int SimulatorBase::PAUSE = -1;
+int SimulatorBase::PAUSE_AT = -1;
+int SimulatorBase::STOP_AT = -1;
+int SimulatorBase::NUM_STEPS_PER_RENDER = -1;
+int SimulatorBase::PARTIO_EXPORT = -1;
+int SimulatorBase::PARTIO_EXPORT_FPS = -1;
+int SimulatorBase::PARTIO_EXPORT_ATTRIBUTES = -1;
+int SimulatorBase::RENDER_WALLS = -1;
+int SimulatorBase::ENUM_WALLS_NONE = -1;
+int SimulatorBase::ENUM_WALLS_PARTICLES_ALL = -1;
+int SimulatorBase::ENUM_WALLS_PARTICLES_NO_WALLS = -1;
+int SimulatorBase::ENUM_WALLS_GEOMETRY_ALL = -1;
+int SimulatorBase::ENUM_WALLS_GEOMETRY_NO_WALLS = -1;
+int SimulatorBase::ENUM_RENDER_NONE = -1;
+int SimulatorBase::ENUM_RENDER_VELOCITY = -1;
+int SimulatorBase::ENUM_RENDER_ANGULAR_VELOCITY = -1;
+int SimulatorBase::ENUM_RENDER_DENSITY = -1;
-DemoBase::DemoBase()
+SimulatorBase::SimulatorBase()
{
Utilities::logger.addSink(unique_ptr<Utilities::ConsoleSink>(new Utilities::ConsoleSink(Utilities::LogLevel::INFO)));
m_numberOfStepsPerRenderUpdate = 4;
- m_colorField = 1;
- setColorMapType(0);
- m_renderMinValue = 0.0;
- m_renderMaxValue = 5.0;
m_sceneFile = "";
m_renderWalls = 4;
m_doPause = true;
@@ -72,17 +63,22 @@ DemoBase::DemoBase()
m_framesPerSecond = 25;
m_nextFrameTime = 0.0;
m_frameCounter = 1;
+ m_colorField.resize(1, "velocity");
+ m_colorMapType.resize(1, 0);
+ m_renderMinValue.resize(1, 0.0);
+ m_renderMaxValue.resize(1, 5.0);
+ m_partioAttributes = "velocity";
#ifdef DL_OUTPUT
m_nextTiming = 1.0;
#endif
}
-DemoBase::~DemoBase()
+SimulatorBase::~SimulatorBase()
{
}
-void DemoBase::initParameters()
+void SimulatorBase::initParameters()
{
ParameterObject::initParameters();
@@ -104,37 +100,10 @@ void DemoBase::initParameters()
setDescription(NUM_STEPS_PER_RENDER, "Number of simulation steps per rendered frame.");
static_cast<NumericParameter<unsigned int>*>(getParameter(NUM_STEPS_PER_RENDER))->setMinValue(1);
- RENDER_MIN_VALUE = createNumericParameter("renderMinValue", "Min. value (shader)", &m_renderMinValue);
- setGroup(RENDER_MIN_VALUE, "Visualization");
- setDescription(RENDER_MIN_VALUE, "Minimal value used for color-coding the color field in the rendering process.");
-
- RENDER_MAX_VALUE = createNumericParameter("renderMaxValue", "Max. value (shader)", &m_renderMaxValue);
- setGroup(RENDER_MAX_VALUE, "Visualization");
- setDescription(RENDER_MAX_VALUE, "Maximal value used for color-coding the color field in the rendering process.");
-
- RENDER_COLOR_FIELD = createEnumParameter("colorField", "Color field", &m_colorField);
- setGroup(RENDER_COLOR_FIELD, "Visualization");
- setDescription(RENDER_COLOR_FIELD, "Choose vector or scalar field for particle coloring.");
- EnumParameter *enumParam = static_cast<EnumParameter*>(getParameter(RENDER_COLOR_FIELD));
- enumParam->addEnumValue("None", ENUM_RENDER_NONE);
- enumParam->addEnumValue("Velocity", ENUM_RENDER_VELOCITY);
- enumParam->addEnumValue("Angular velocity (micropolar model)", ENUM_RENDER_ANGULAR_VELOCITY);
- enumParam->addEnumValue("Density", ENUM_RENDER_DENSITY);
-
- ParameterBase::GetFunc<int> getColorMapTypeFct = std::bind(&DemoBase::getColorMapType, this);
- ParameterBase::SetFunc<int> setColorMapTypeFct = std::bind(&DemoBase::setColorMapType, this, std::placeholders::_1);
- RENDER_COLOR_MAP_TYPE = createEnumParameter("colorMapType", "Color map", getColorMapTypeFct, setColorMapTypeFct);
- setGroup(RENDER_COLOR_MAP_TYPE, "Visualization");
- setDescription(RENDER_COLOR_MAP_TYPE, "Choose a color map.");
- enumParam = static_cast<EnumParameter*>(getParameter(RENDER_COLOR_MAP_TYPE));
- enumParam->addEnumValue("None", ENUM_COLORMAP_NONE);
- enumParam->addEnumValue("Jet", ENUM_COLORMAP_JET);
- enumParam->addEnumValue("Plasma", ENUM_COLORMAP_PLASMA);
-
RENDER_WALLS = createEnumParameter("renderWalls", "Render walls", &m_renderWalls);
setGroup(RENDER_WALLS, "Visualization");
setDescription(RENDER_WALLS, "Make walls visible/invisible.");
- enumParam = static_cast<EnumParameter*>(getParameter(RENDER_WALLS));
+ EnumParameter *enumParam = static_cast<EnumParameter*>(getParameter(RENDER_WALLS));
enumParam->addEnumValue("None", ENUM_WALLS_NONE);
enumParam->addEnumValue("Particles (all)", ENUM_WALLS_PARTICLES_ALL);
enumParam->addEnumValue("Particles (no walls)", ENUM_WALLS_PARTICLES_NO_WALLS);
@@ -148,9 +117,14 @@ void DemoBase::initParameters()
PARTIO_EXPORT_FPS = createNumericParameter("partioFPS", "Export FPS", &m_framesPerSecond);
setGroup(PARTIO_EXPORT_FPS, "Export");
setDescription(PARTIO_EXPORT_FPS, "Frame rate of partio export.");
+
+ PARTIO_EXPORT_ATTRIBUTES = createStringParameter("partioAttributes", "Export attributes", &m_partioAttributes);
+ getParameter(PARTIO_EXPORT_ATTRIBUTES)->setReadOnly(true);
+ setGroup(PARTIO_EXPORT_ATTRIBUTES, "Export");
+ setDescription(PARTIO_EXPORT_ATTRIBUTES, "Attributes that are exported in the partio files (except id and position).");
}
-void DemoBase::init(int argc, char **argv, const char *demoName)
+void SimulatorBase::init(int argc, char **argv, const char *simName)
{
initParameters();
m_exePath = FileSystem::getProgramPath();
@@ -199,18 +173,22 @@ void DemoBase::init(int argc, char **argv, const char *demoName)
return;
// OpenGL
- MiniGL::init(argc, argv, 1280, 960, 0, 0, demoName);
+ MiniGL::init(argc, argv, 1280, 960, 0, 0, simName);
MiniGL::initLights();
MiniGL::getOpenGLVersion(m_context_major_version, m_context_minor_version);
- MiniGL::setViewport(40.0, 0.1f, 500.0, Vector3r(0.0, 3.0, 8.0), Vector3r(0.0, 0.0, 0.0));
+ const bool sim2D = getScene().sim2D;
+ if (sim2D)
+ MiniGL::setViewport(40.0, 0.1f, 500.0, Vector3r(0.0, 0.0, 8.0), Vector3r(0.0, 0.0, 0.0));
+ else
+ MiniGL::setViewport(40.0, 0.1f, 500.0, Vector3r(0.0, 3.0, 8.0), Vector3r(0.0, 0.0, 0.0));
MiniGL::setSelectionFunc(selection, this);
- MiniGL::addKeyFunc('i', std::bind(&DemoBase::particleInfo, this));
+ MiniGL::addKeyFunc('i', std::bind(&SimulatorBase::particleInfo, this));
if (MiniGL::checkOpenGLVersion(3, 3))
initShaders();
}
-void DemoBase::cleanup()
+void SimulatorBase::cleanup()
{
for (unsigned int i = 0; i < m_scene.boundaryModels.size(); i++)
delete m_scene.boundaryModels[i];
@@ -229,7 +207,7 @@ void DemoBase::cleanup()
m_scene.emitters.clear();
}
-void DemoBase::initShaders()
+void SimulatorBase::initShaders()
{
string vertFile = getExePath() + "/resources/shaders/vs_points_vector.glsl";
string fragFile = getExePath() + "/resources/shaders/fs_points.glsl";
@@ -305,7 +283,7 @@ void DemoBase::initShaders()
}
-void DemoBase::meshShaderBegin(const float *col)
+void SimulatorBase::meshShaderBegin(const float *col)
{
m_meshShader.begin();
glUniform1f(m_meshShader.getUniform("shininess"), 5.0f);
@@ -320,12 +298,12 @@ void DemoBase::meshShaderBegin(const float *col)
glUniform3fv(m_meshShader.getUniform("surface_color"), 1, col);
}
-void DemoBase::meshShaderEnd()
+void SimulatorBase::meshShaderEnd()
{
m_meshShader.end();
}
-void DemoBase::pointShaderBegin(Shader *shader, const float *col, const bool useTexture)
+void SimulatorBase::pointShaderBegin(Shader *shader, const float *col, const Real minVal, const Real maxVal, const bool useTexture, float const* color_map)
{
shader->begin();
@@ -334,16 +312,16 @@ void DemoBase::pointShaderBegin(Shader *shader, const float *col, const bool use
const Real radius = Simulation::getCurrent()->getValue<Real>(Simulation::PARTICLE_RADIUS);
glUniform1f(shader->getUniform("viewport_width"), (float)viewport[2]);
glUniform1f(shader->getUniform("radius"), (float)radius);
- glUniform1f(shader->getUniform("min_scalar"), (GLfloat)m_renderMinValue);
- glUniform1f(shader->getUniform("max_scalar"), (GLfloat)m_renderMaxValue);
+ glUniform1f(shader->getUniform("min_scalar"), (GLfloat)minVal);
+ glUniform1f(shader->getUniform("max_scalar"), (GLfloat)maxVal);
glUniform3fv(shader->getUniform("color"), 1, col);
if (useTexture)
{
glActiveTexture(GL_TEXTURE0);
glBindTexture(GL_TEXTURE_1D, m_textureMap);
- glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, m_colorMapLength, 0, GL_RGB, GL_FLOAT,
- reinterpret_cast<float const*>(m_colorMapBuffer));
+ glTexImage1D(GL_TEXTURE_1D, 0, GL_RGB, 256u, 0, GL_RGB, GL_FLOAT, color_map);
+
glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MAG_FILTER, GL_NEAREST);
glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_MIN_FILTER, GL_NEAREST);
glTexParameteri(GL_TEXTURE_1D, GL_TEXTURE_WRAP_S, GL_CLAMP_TO_EDGE);
@@ -367,13 +345,13 @@ void DemoBase::pointShaderBegin(Shader *shader, const float *col, const bool use
glPointParameterf(GL_POINT_SPRITE_COORD_ORIGIN, GL_LOWER_LEFT);
}
-void DemoBase::pointShaderEnd(Shader *shader, const bool useTexture)
+void SimulatorBase::pointShaderEnd(Shader *shader, const bool useTexture)
{
glBindTexture(GL_TEXTURE_1D, 0);
shader->end();
}
-void DemoBase::readParameters()
+void SimulatorBase::readParameters()
{
m_sceneLoader->readParameterObject("Configuration", this);
m_sceneLoader->readParameterObject("Configuration", Simulation::getCurrent());
@@ -389,11 +367,18 @@ void DemoBase::readParameters()
m_sceneLoader->readParameterObject(key, (ParameterObject*) model->getSurfaceTensionBase());
m_sceneLoader->readParameterObject(key, (ParameterObject*) model->getViscosityBase());
m_sceneLoader->readParameterObject(key, (ParameterObject*) model->getVorticityBase());
+ m_sceneLoader->readParameterObject(key, (ParameterObject*) model->getElasticityBase());
+
+ SceneLoader::ColoringData colorData = m_sceneLoader->readColoringInfo(key);
+ setColorField(i, colorData.colorField);
+ setColorMapType(i, colorData.colorMapType);
+ setRenderMinValue(i, colorData.minVal);
+ setRenderMaxValue(i, colorData.maxVal);
}
}
-void DemoBase::buildModel()
+void SimulatorBase::buildModel()
{
TimeManager::getCurrent()->setTimeStepSize(m_scene.timeStepSize);
@@ -403,7 +388,8 @@ void DemoBase::buildModel()
Simulation *sim = Simulation::getCurrent();
- sim->getTimeStep()->resize();
+ if (sim->getTimeStep())
+ sim->getTimeStep()->resize();
if (!sim->is2DSimulation())
{
@@ -418,7 +404,7 @@ void DemoBase::buildModel()
}
-void DemoBase::initFluidData()
+void SimulatorBase::initFluidData()
{
LOG_INFO << "Initialize fluid particles";
@@ -430,7 +416,7 @@ void DemoBase::initFluidData()
std::map<std::string, unsigned int> fluidIDs;
unsigned int index = 0;
for (unsigned int i = 0; i < m_scene.fluidBlocks.size(); i++)
- {
+ {
if (fluidIDs.find(m_scene.fluidBlocks[i]->id) == fluidIDs.end())
fluidIDs[m_scene.fluidBlocks[i]->id] = index++;
}
@@ -482,10 +468,15 @@ void DemoBase::initFluidData()
nParticles += (unsigned int)fluidParticles[index].size();
}
+ m_colorField.resize(sim->numberOfFluidModels(), "velocity");
+ m_colorMapType.resize(sim->numberOfFluidModels(), 0);
+ m_renderMinValue.resize(sim->numberOfFluidModels(), 0.0);
+ m_renderMaxValue.resize(sim->numberOfFluidModels(), 5.0);
+
LOG_INFO << "Number of fluid particles: " << nParticles;
}
-void DemoBase::createEmitters()
+void SimulatorBase::createEmitters()
{
Simulation *sim = Simulation::getCurrent();
@@ -531,7 +522,7 @@ void DemoBase::createEmitters()
}
-void DemoBase::createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs, std::vector<std::vector<Vector3r>> &fluidParticles, std::vector<std::vector<Vector3r>> &fluidVelocities)
+void SimulatorBase::createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs, std::vector<std::vector<Vector3r>> &fluidParticles, std::vector<std::vector<Vector3r>> &fluidVelocities)
{
for (unsigned int i = 0; i < m_scene.fluidBlocks.size(); i++)
{
@@ -611,9 +602,11 @@ void DemoBase::createFluidBlocks(std::map<std::string, unsigned int> &fluidIDs,
}
}
-void DemoBase::renderFluid(FluidModel *model, float *fluidColor)
+void SimulatorBase::renderFluid(const unsigned int fluidModelIndex, float *fluidColor)
{
// Draw simulation model
+ Simulation *sim = Simulation::getCurrent();
+ FluidModel *model = sim->getFluidModel(fluidModelIndex);
const unsigned int nParticles = model->numActiveParticles();
if (nParticles == 0)
return;
@@ -625,8 +618,7 @@ void DemoBase::renderFluid(FluidModel *model, float *fluidColor)
glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, speccolor);
glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 100.0);
glColor3fv(surfaceColor);
-
- Simulation *sim = Simulation::getCurrent();
+
const Real supportRadius = sim->getSupportRadius();
Real vmax = static_cast<Real>(0.4*2.0)*supportRadius / TimeManager::getCurrent()->getTimeStepSize();
Real vmin = 0.0;
@@ -635,56 +627,60 @@ void DemoBase::renderFluid(FluidModel *model, float *fluidColor)
{
Shader *shader_vector = &m_shader_vector_map;
Shader *shader_scalar = &m_shader_scalar_map;
- if (m_colorMapType == ENUM_COLORMAP_NONE)
+ float const *color_map = nullptr;
+ if (m_colorMapType[fluidModelIndex] == 1)
+ color_map = reinterpret_cast<float const*>(colormap_jet);
+ else if (m_colorMapType[fluidModelIndex] == 2)
+ color_map = reinterpret_cast<float const*>(colormap_plasma);
+
+ if (m_colorMapType[fluidModelIndex] == 0)
{
shader_vector = &m_shader_vector;
shader_scalar = &m_shader_scalar;
}
- if ((m_colorField == ENUM_RENDER_VELOCITY) || (m_colorField == ENUM_RENDER_ANGULAR_VELOCITY))
- pointShaderBegin(shader_vector, &fluidColor[0], true);
- else if (m_colorField == ENUM_RENDER_DENSITY)
- pointShaderBegin(shader_scalar, &fluidColor[0], true);
- else
- pointShaderBegin(shader_scalar, &fluidColor[0], false);
+ const FieldDescription *field = nullptr;
+ const std::string &colorFieldName = m_colorField[fluidModelIndex];
+ field = &model->getField(colorFieldName);
+
+
+ if (field == nullptr)
+ pointShaderBegin(shader_scalar, &fluidColor[0], m_renderMinValue[fluidModelIndex], m_renderMaxValue[fluidModelIndex], false);
+ else if (field->type == FieldType::Vector3)
+ pointShaderBegin(shader_vector, &fluidColor[0], m_renderMinValue[fluidModelIndex], m_renderMaxValue[fluidModelIndex], true, color_map);
+ else if (field->type == FieldType::Scalar)
+ pointShaderBegin(shader_scalar, &fluidColor[0], m_renderMinValue[fluidModelIndex], m_renderMaxValue[fluidModelIndex], true, color_map);
if (model->numActiveParticles() > 0)
{
glEnableVertexAttribArray(0);
glVertexAttribPointer(0, 3, GL_REAL, GL_FALSE, 0, &model->getPosition(0));
-
- if (m_colorField == ENUM_RENDER_VELOCITY)
- {
- glEnableVertexAttribArray(1);
- glVertexAttribPointer(1, 3, GL_REAL, GL_FALSE, 0, &model->getVelocity(0));
- }
- else if ((m_colorField == ENUM_RENDER_ANGULAR_VELOCITY) && ((VorticityMethods)model->getVorticityMethod() == VorticityMethods::Micropolar))
- {
- glEnableVertexAttribArray(1);
- float fluidColor2[4] = { 0.3f, 0.9f, 0.5f, 1.0f };
- glUniform3fv(shader_vector->getUniform("color"), 1, fluidColor2);
- glVertexAttribPointer(1, 3, GL_REAL, GL_FALSE, 0, &((MicropolarModel_Bender2017*)model->getVorticityBase())->getAngularVelocity(0)[0]);
- }
- else if (m_colorField == ENUM_RENDER_DENSITY)
+
+ if (field != nullptr)
{
- glEnableVertexAttribArray(1);
- float fluidColor2[4] = { 0.9f, 0.3f, 0.3f, 1.0f };
- glUniform3fv(shader_scalar->getUniform("color"), 1, fluidColor2);
- glVertexAttribPointer(1, 1, GL_REAL, GL_FALSE, 0, &model->getDensity(0));
- }
+ if (field->type == FieldType::Vector3)
+ {
+ glEnableVertexAttribArray(1);
+ glVertexAttribPointer(1, 3, GL_REAL, GL_FALSE, 0, field->getFct(0));
+ }
+ else if (field->type == FieldType::Scalar)
+ {
+ glEnableVertexAttribArray(1);
+ glVertexAttribPointer(1, 1, GL_REAL, GL_FALSE, 0, field->getFct(0));
+ }
+ }
glDrawArrays(GL_POINTS, 0, model->numActiveParticles());
glDisableVertexAttribArray(0);
glDisableVertexAttribArray(1);
}
- if ((m_colorField == ENUM_RENDER_VELOCITY) || (m_colorField == ENUM_RENDER_ANGULAR_VELOCITY))
+ if (field == nullptr)
+ pointShaderEnd(shader_scalar, false);
+ else if (field->type == FieldType::Vector3)
pointShaderEnd(shader_vector, true);
- else if (m_colorField == ENUM_RENDER_DENSITY)
+ else if (field->type == FieldType::Scalar)
pointShaderEnd(shader_scalar, true);
- else
- pointShaderEnd(shader_scalar, false);
-
}
else
{
@@ -711,7 +707,7 @@ void DemoBase::renderFluid(FluidModel *model, float *fluidColor)
const unsigned int fluidIndex = model->getPointSetIndex();
if (MiniGL::checkOpenGLVersion(3, 3))
{
- pointShaderBegin(&m_shader_scalar, &red[0]);
+ pointShaderBegin(&m_shader_scalar, &red[0], m_renderMinValue[fluidModelIndex], m_renderMaxValue[fluidModelIndex]);
if ((getSelectedParticles().size() > 0) && ((getSelectedParticles()[fluidIndex].size() > 0)))
{
const Real radius = sim->getValue<Real>(Simulation::PARTICLE_RADIUS);
@@ -745,9 +741,9 @@ void DemoBase::renderFluid(FluidModel *model, float *fluidColor)
}
-void DemoBase::mouseMove(int x, int y, void *clientData)
+void SimulatorBase::mouseMove(int x, int y, void *clientData)
{
- DemoBase *base = (DemoBase*)clientData;
+ SimulatorBase *base = (SimulatorBase*)clientData;
Simulation *sim = Simulation::getCurrent();
Vector3r mousePos;
@@ -768,9 +764,9 @@ void DemoBase::mouseMove(int x, int y, void *clientData)
base->m_oldMousePos = mousePos;
}
-void DemoBase::selection(const Eigen::Vector2i &start, const Eigen::Vector2i &end, void *clientData)
+void SimulatorBase::selection(const Eigen::Vector2i &start, const Eigen::Vector2i &end, void *clientData)
{
- DemoBase *base = (DemoBase*)clientData;
+ SimulatorBase *base = (SimulatorBase*)clientData;
Simulation *sim = Simulation::getCurrent();
base->m_selectedParticles.resize(sim->numberOfFluidModels());
bool selected = false;
@@ -803,24 +799,59 @@ void DemoBase::selection(const Eigen::Vector2i &start, const Eigen::Vector2i &en
}
}
-void DemoBase::particleInfo()
+void SimulatorBase::particleInfo()
{
Simulation *sim = Simulation::getCurrent();
+ const int maxWidth = 25;
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i++)
{
FluidModel *model = sim->getFluidModel(i);
+ if (m_selectedParticles[i].size() > 0)
+ {
+ LOG_INFO << "---------------------------------------------------------------------------";
+ LOG_INFO << model->getId();
+ LOG_INFO << "---------------------------------------------------------------------------";
+ }
for (unsigned int j = 0; j < m_selectedParticles[i].size(); j++)
{
unsigned int index = m_selectedParticles[i][j];
- LOG_INFO << index;
- LOG_INFO << "x: " << model->getPosition(index).transpose();
- LOG_INFO << "v: " << model->getVelocity(index).transpose();
- LOG_INFO << "density: " << model->getDensity(index);
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << "Index:" << index;
+ for (unsigned int k = 0; k < model->numberOfFields(); k++)
+ {
+ const FieldDescription &field = model->getField(k);
+ if (field.type == Scalar)
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name + ":" << *field.getFct(index);
+ else if (field.type == Vector3)
+ {
+ Eigen::Map<Vector3r> vec(field.getFct(index));
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name + ":" << vec.transpose();
+ }
+ else if (field.type == Vector6)
+ {
+ Eigen::Map<Vector6r> vec(field.getFct(index));
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name + ":" << vec.transpose();
+ }
+ else if (field.type == Matrix3)
+ {
+ Eigen::Map<Matrix3r> mat(field.getFct(index));
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name + ":" << mat.row(0);
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(1);
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(2);
+ }
+ else if (field.type == Matrix6)
+ {
+ Eigen::Map<Matrix6r> mat(field.getFct(index));
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << field.name + ":" << mat.row(0);
+ for (unsigned int k = 1; k < 6; k++)
+ LOG_INFO << std::left << std::setw(maxWidth) << std::setfill(' ') << " " << mat.row(k);
+ }
+ }
+ LOG_INFO << "---------------------------------------------------------------------------\n";
}
}
}
-void DemoBase::partioExport()
+void SimulatorBase::partioExport()
{
std::string exportPath = FileSystem::normalizePath(m_outputPath + "/partio");
FileSystem::makeDirs(exportPath);
@@ -833,11 +864,107 @@ void DemoBase::partioExport()
fileName = fileName + "_" + model->getId() + "_" + std::to_string(m_frameCounter) + ".bgeo";
std::string exportFileName = FileSystem::normalizePath(exportPath + "/" + fileName);
- PartioReaderWriter::writeParticles(exportFileName, model->numActiveParticles(), &model->getPosition(0), &model->getVelocity(0), 0.0);
+ writeParticles(exportFileName, model);
}
}
-void DemoBase::step()
+
+void SimulatorBase::writeParticles(const std::string &fileName, FluidModel *model)
+{
+ Partio::ParticlesDataMutable& particleData = *Partio::create();
+ Partio::ParticleAttribute posAttr = particleData.addAttribute("position", Partio::VECTOR, 3);
+ Partio::ParticleAttribute idAttr = particleData.addAttribute("id", Partio::INT, 1);
+
+ // add attributes
+ std::vector<std::string> attributes;
+ StringTools::tokenize(m_partioAttributes, attributes, ";");
+
+ std::map<unsigned int, int> attrMap;
+ std::map<unsigned int, Partio::ParticleAttribute> partioAttrMap;
+ for (unsigned int i = 0; i < attributes.size(); i++)
+ {
+ // position is exported anyway
+ if (attributes[i] == "position")
+ {
+ attrMap[i] = -1;
+ continue;
+ }
+
+ bool found = false;
+ for (unsigned int j = 0; j < model->numberOfFields(); j++)
+ {
+ const FieldDescription &field = model->getField(j);
+ if (field.name == attributes[i])
+ {
+ found = true;
+ if (field.type == Scalar)
+ {
+ attrMap[i] = j;
+ partioAttrMap[i] = particleData.addAttribute(attributes[i].c_str(), Partio::FLOAT, 1);
+ }
+ else if (field.type == Vector3)
+ {
+ attrMap[i] = j;
+ partioAttrMap[i] = particleData.addAttribute(attributes[i].c_str(), Partio::VECTOR, 3);
+ }
+ else
+ {
+ attrMap[i] = -1;
+ LOG_WARN << "Only scalar and vector fields are currently supported by the partio exporter.";
+ }
+ break;
+ }
+ }
+ if (!found)
+ {
+ attrMap[i] = -1;
+ LOG_WARN << "Unknown field cannot be exported in partio file: " << attributes[i];
+ }
+ }
+
+ const unsigned int numParticles = model->numActiveParticles();
+
+ for (unsigned int i = 0; i < numParticles; i++)
+ {
+ Partio::ParticleIndex index = particleData.addParticle();
+ float* pos = particleData.dataWrite<float>(posAttr, index);
+ int* id = particleData.dataWrite<int>(idAttr, index);
+
+ const Vector3r &x = model->getPosition(i);
+ pos[0] = (float)x[0];
+ pos[1] = (float)x[1];
+ pos[2] = (float)x[2];
+
+ id[0] = i;
+
+ for (unsigned int j = 0; j < attributes.size(); j++)
+ {
+ const int fieldIndex = attrMap[j];
+ if (fieldIndex != -1)
+ {
+ const FieldDescription &field = model->getField(fieldIndex);
+ if (field.type == FieldType::Scalar)
+ {
+ float* val = particleData.dataWrite<float>(partioAttrMap[j], index);
+ *val = (float) *field.getFct(i);
+ }
+ else if (field.type == FieldType::Vector3)
+ {
+ float* val = particleData.dataWrite<float>(partioAttrMap[j], index);
+ Eigen::Map<Vector3r> vec(field.getFct(i));
+ val[0] = (float)vec[0];
+ val[1] = (float)vec[1];
+ val[2] = (float)vec[2];
+ }
+ }
+ }
+ }
+
+ Partio::write(fileName.c_str(), particleData, true);
+ particleData.release();
+}
+
+void SimulatorBase::step()
{
if (m_enablePartioExport)
{
@@ -862,7 +989,7 @@ void DemoBase::step()
#endif
}
-void DemoBase::reset()
+void SimulatorBase::reset()
{
//TimeManager::getCurrent()->setTimeStepSize(m_scene.timeStepSize);
m_nextFrameTime = 0.0;
@@ -872,17 +999,3 @@ void DemoBase::reset()
#endif
}
-void DemoBase::setColorMapType(const int v)
-{
- m_colorMapType = v;
- if (m_colorMapType == ENUM_COLORMAP_JET)
- {
- m_colorMapBuffer = colormap_jet[0];
- m_colorMapLength = 256u;
- }
- else if (m_colorMapType == ENUM_COLORMAP_PLASMA)
- {
- m_colorMapBuffer = colormap_plasma[0];
- m_colorMapLength = 256u;
- }
-}
diff --git a/Demos/Common/DemoBase.h b/Simulators/Common/SimulatorBase.h
similarity index 68%
rename from Demos/Common/DemoBase.h
rename to Simulators/Common/SimulatorBase.h
index 4d6ea9b..cb0dd4a 100644
--- a/Demos/Common/DemoBase.h
+++ b/Simulators/Common/SimulatorBase.h
@@ -1,5 +1,5 @@
-#ifndef __DemoBase_h__
-#define __DemoBase_h__
+#ifndef __SimulatorBase_h__
+#define __SimulatorBase_h__
#include "SPlisHSPlasH/Common.h"
#include "SPlisHSPlasH/Utilities/SceneLoader.h"
@@ -11,7 +11,7 @@
namespace SPH
{
- class DemoBase : public GenParam::ParameterObject
+ class SimulatorBase : public GenParam::ParameterObject
{
public:
struct SimulationMethod
@@ -38,20 +38,21 @@ namespace SPH
Shader m_meshShader;
GLuint m_textureMap;
int m_renderWalls;
- int m_colorField;
bool m_doPause;
Real m_pauseAt;
Real m_stopAt;
bool m_enablePartioExport;
unsigned int m_framesPerSecond;
- Real m_renderMaxValue;
- Real m_renderMinValue;
+ std::string m_partioAttributes;
Vector3r m_oldMousePos;
std::vector<std::vector<unsigned int>> m_selectedParticles;
std::unique_ptr<Utilities::SceneLoader> m_sceneLoader;
Real m_nextFrameTime;
unsigned int m_frameCounter;
- int m_colorMapType;
+ std::vector<std::string> m_colorField;
+ std::vector<int> m_colorMapType;
+ std::vector<Real> m_renderMaxValue;
+ std::vector<Real> m_renderMinValue;
float const* m_colorMapBuffer;
unsigned int m_colorMapLength;
#ifdef DL_OUTPUT
@@ -76,12 +77,9 @@ namespace SPH
static int NUM_STEPS_PER_RENDER;
static int PARTIO_EXPORT;
static int PARTIO_EXPORT_FPS;
- static int RENDER_MIN_VALUE;
- static int RENDER_MAX_VALUE;
+ static int PARTIO_EXPORT_ATTRIBUTES;
static int RENDER_WALLS;
- static int RENDER_COLOR_FIELD;
- static int RENDER_COLOR_MAP_TYPE;
-
+
static int ENUM_WALLS_NONE;
static int ENUM_WALLS_PARTICLES_ALL;
static int ENUM_WALLS_PARTICLES_NO_WALLS;
@@ -93,21 +91,18 @@ namespace SPH
static int ENUM_RENDER_ANGULAR_VELOCITY;
static int ENUM_RENDER_DENSITY;
- static int ENUM_COLORMAP_NONE;
- static int ENUM_COLORMAP_JET;
- static int ENUM_COLORMAP_PLASMA;
-
- DemoBase();
- virtual ~DemoBase();
+ SimulatorBase();
+ virtual ~SimulatorBase();
- void init(int argc, char **argv, const char *demoName);
+ void init(int argc, char **argv, const char *simName);
void buildModel();
void cleanup();
- void renderFluid(FluidModel *model, float *fluidColor);
+ void renderFluid(const unsigned int fluidModelIndex, float *fluidColor);
void readParameters();
void partioExport();
+ void writeParticles(const std::string &fileName, FluidModel *model);
void step();
void reset();
@@ -124,7 +119,7 @@ namespace SPH
Shader& getMeshShader() { return m_meshShader; }
void meshShaderBegin(const float *col);
void meshShaderEnd();
- void pointShaderBegin(Shader *shader, const float *col, const bool useTexture = false);
+ void pointShaderBegin(Shader *shader, const float *col, const Real minVal, const Real maxVal, const bool useTexture = false, float const* color_map = nullptr);
void pointShaderEnd(Shader *shader, const bool useTexture = false);
Utilities::SceneLoader::Scene& getScene() { return m_scene; }
@@ -132,8 +127,15 @@ namespace SPH
bool getUseParticleCaching() const { return m_useParticleCaching; }
void setUseParticleCaching(bool val) { m_useParticleCaching = val; }
- int getColorMapType() const { return m_colorMapType; }
- void setColorMapType(const int v);
+ 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; }
+
+ int getColorMapType(const unsigned int fluidModelIndex) const { return m_colorMapType[fluidModelIndex]; }
+ void setColorMapType(const unsigned int fluidModelIndex, int val) { m_colorMapType[fluidModelIndex] = val; }
+ Real getRenderMaxValue(const unsigned int fluidModelIndex) const { return m_renderMaxValue[fluidModelIndex]; }
+ void setRenderMaxValue(const unsigned int fluidModelIndex, Real val) { m_renderMaxValue[fluidModelIndex] = val; }
+ Real getRenderMinValue(const unsigned int fluidModelIndex) const { return m_renderMinValue[fluidModelIndex]; }
+ void setRenderMinValue(const unsigned int fluidModelIndex, Real val) { m_renderMinValue[fluidModelIndex] = val; }
};
}
diff --git a/Demos/Common/TweakBarParameters.cpp b/Simulators/Common/TweakBarParameters.cpp
similarity index 89%
rename from Demos/Common/TweakBarParameters.cpp
rename to Simulators/Common/TweakBarParameters.cpp
index fae8234..c51607d 100644
--- a/Demos/Common/TweakBarParameters.cpp
+++ b/Simulators/Common/TweakBarParameters.cpp
@@ -130,25 +130,11 @@ void TweakBarParameters::createParameterObjectGUI(ParameterObject *paramObj)
TwAddVarCB(MiniGL::getTweakBar(), varName, TW_TYPE_DIR3R, setParameterValue, getParameterValue, m_params[m_params.size() - 1].get(), str);
}
}
-
-// else if ((desc->type == IBDS::IBDS_STRING_PARAMETER) || (desc->type == IBDS::IBDS_FILEPATH_PARAMETER))
-// {
-// //const bool &defaultValue = paramObj->getParameterDefaultValue<bool>(j);
-// char str[400];
-// char tmp.c_str()[20];
-// if (desc->readOnly)
-// sprintf(tmp.c_str(), "%s", "readonly=true");
-// else
-// sprintf(tmp.c_str(), "%s", "readonly=false");
-//
-// // Defining new unique var name
-// char varName[500];
-// sprintf(varName, "%s::%s", objName.c_str(), desc->fullName.c_str());
-//
-// sprintf(str, " label='%s' help='%s' group='%s' %s ", desc->label.c_str(), desc->description.c_str(), objName.c_str(), tmp.c_str());
-// TwAddVarCB(MiniGL::getTwBar(), varName, TW_TYPE_CDSTRING, setParameterValue, getParameterValue, &m_params[m_params.size() - 1], str);
-// }
-//
+ else if (paramBase->getType() == ParameterBase::STRING)
+ {
+ sprintf(str, " label='%s' help='%s' group='%s' %s ", paramBase->getLabel().c_str(), paramBase->getDescription().c_str(), paramBase->getGroup().c_str(), tmp.c_str());
+ TwAddVarCB(MiniGL::getTweakBar(), varName, TW_TYPE_STDSTRING, setParameterValue, getParameterValue, m_params[m_params.size() - 1].get(), str);
+ }
}
}
@@ -232,11 +218,11 @@ void TW_CALL TweakBarParameters::setParameterValue(const void *value, void *clie
static_cast<VectorParameter<Real>*>(paramBase)->setValue(val);
}
}
-// else if ((desc->type == IBDS::IBDS_STRING_PARAMETER) || (desc->type == IBDS::IBDS_FILEPATH_PARAMETER))
-// {
-// const std::string val = std::string(*(const char **)(value));
-// paramObj->setParameter<std::string>(desc->id, val);
-// }
+ else if (paramBase->getType() == ParameterBase::STRING)
+ {
+ const std::string *srcPtr = static_cast<const std::string *>(value);
+ static_cast<StringParameter*>(paramBase)->setValue(*srcPtr);
+ }
}
void TW_CALL TweakBarParameters::getParameterValue(void *value, void *clientData)
@@ -300,10 +286,11 @@ void TW_CALL TweakBarParameters::getParameterValue(void *value, void *clientData
((Real*)value)[2] = val[2];
}
}
-// else if ((desc->type == IBDS::IBDS_STRING_PARAMETER) || (desc->type == IBDS::IBDS_FILEPATH_PARAMETER))
-// {
-// const std::string &val = paramObj->getParameter<std::string>(desc->id);
-// *(const char **)(value) = val.c_str();
-// }
+ else if (paramBase->getType() == ParameterBase::STRING)
+ {
+ const std::string &val = static_cast<StringParameter*>(paramBase)->getValue();
+ std::string *destPtr = static_cast<std::string *>(value);
+ TwCopyStdStringToLibrary(*destPtr, val);
+ }
}
diff --git a/Demos/Common/TweakBarParameters.h b/Simulators/Common/TweakBarParameters.h
similarity index 100%
rename from Demos/Common/TweakBarParameters.h
rename to Simulators/Common/TweakBarParameters.h
diff --git a/Demos/DynamicBoundaryDemo/CMakeLists.txt b/Simulators/DynamicBoundarySimulator/CMakeLists.txt
similarity index 75%
rename from Demos/DynamicBoundaryDemo/CMakeLists.txt
rename to Simulators/DynamicBoundarySimulator/CMakeLists.txt
index 69d8e65..cd7e6a4 100644
--- a/Demos/DynamicBoundaryDemo/CMakeLists.txt
+++ b/Simulators/DynamicBoundarySimulator/CMakeLists.txt
@@ -64,13 +64,13 @@ set(SIMULATION_LINK_LIBRARIES ${SIMULATION_LINK_LIBRARIES}
set(SIMULATION_DEPENDENCIES Ext_PBD ${SIMULATION_DEPENDENCIES})
link_directories(${PROJECT_PATH}/extern/install/PositionBasedDynamics/lib)
-add_executable(DynamicBoundaryDemo
+add_executable(DynamicBoundarySimulator
main.cpp
- ${PROJECT_PATH}/Demos/Common/DemoBase.cpp
- ${PROJECT_PATH}/Demos/Common/DemoBase.h
- ${PROJECT_PATH}/Demos/Common/TweakBarParameters.cpp
- ${PROJECT_PATH}/Demos/Common/TweakBarParameters.h
+ ${PROJECT_PATH}/Simulators/Common/SimulatorBase.cpp
+ ${PROJECT_PATH}/Simulators/Common/SimulatorBase.h
+ ${PROJECT_PATH}/Simulators/Common/TweakBarParameters.cpp
+ ${PROJECT_PATH}/Simulators/Common/TweakBarParameters.h
${VIS_FILES}
@@ -87,23 +87,23 @@ add_definitions(-DTW_NO_LIB_PRAGMA -DTW_STATIC)
include_directories(${PROJECT_PATH}/extern/freeglut/include)
include_directories(${PROJECT_PATH}/extern/glew/include)
-set_target_properties(DynamicBoundaryDemo PROPERTIES DEBUG_POSTFIX ${CMAKE_DEBUG_POSTFIX})
-set_target_properties(DynamicBoundaryDemo PROPERTIES RELWITHDEBINFO_POSTFIX ${CMAKE_RELWITHDEBINFO_POSTFIX})
-set_target_properties(DynamicBoundaryDemo PROPERTIES MINSIZEREL_POSTFIX ${CMAKE_MINSIZEREL_POSTFIX})
-add_dependencies(DynamicBoundaryDemo ${SIMULATION_DEPENDENCIES})
-target_link_libraries(DynamicBoundaryDemo ${SIMULATION_LINK_LIBRARIES})
+set_target_properties(DynamicBoundarySimulator PROPERTIES DEBUG_POSTFIX ${CMAKE_DEBUG_POSTFIX})
+set_target_properties(DynamicBoundarySimulator PROPERTIES RELWITHDEBINFO_POSTFIX ${CMAKE_RELWITHDEBINFO_POSTFIX})
+set_target_properties(DynamicBoundarySimulator PROPERTIES MINSIZEREL_POSTFIX ${CMAKE_MINSIZEREL_POSTFIX})
+add_dependencies(DynamicBoundarySimulator ${SIMULATION_DEPENDENCIES})
+target_link_libraries(DynamicBoundarySimulator ${SIMULATION_LINK_LIBRARIES})
VIS_SOURCE_GROUPS()
source_group("Header Files\\PBD" FILES ${PBDWRAPPER_HEADER_FILES})
source_group("Source Files\\PBD" FILES ${PBDWRAPPER_SOURCE_FILES})
-set_target_properties(DynamicBoundaryDemo PROPERTIES FOLDER "Demos")
+set_target_properties(DynamicBoundarySimulator PROPERTIES FOLDER "Simulators")
-add_custom_command(TARGET DynamicBoundaryDemo PRE_BUILD
+add_custom_command(TARGET DynamicBoundarySimulator PRE_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_directory
- ${PROJECT_PATH}/data/shaders $<TARGET_FILE_DIR:DynamicBoundaryDemo>/resources/shaders)
+ ${PROJECT_PATH}/data/shaders $<TARGET_FILE_DIR:DynamicBoundarySimulator>/resources/shaders)
-add_custom_command(TARGET DynamicBoundaryDemo PRE_BUILD
+add_custom_command(TARGET DynamicBoundarySimulator PRE_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_directory
- ${PROJECT_PATH}/extern/install/PositionBasedDynamics/include/data/shaders $<TARGET_FILE_DIR:DynamicBoundaryDemo>/resources/pbd_shaders)
+ ${PROJECT_PATH}/extern/install/PositionBasedDynamics/include/data/shaders $<TARGET_FILE_DIR:DynamicBoundarySimulator>/resources/pbd_shaders)
diff --git a/Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDRigidBody.h b/Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDRigidBody.h
similarity index 100%
rename from Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDRigidBody.h
rename to Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDRigidBody.h
diff --git a/Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDWrapper.cpp b/Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDWrapper.cpp
similarity index 99%
rename from Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDWrapper.cpp
rename to Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDWrapper.cpp
index 34843b8..9241792 100644
--- a/Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDWrapper.cpp
+++ b/Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDWrapper.cpp
@@ -10,10 +10,10 @@
#include "Utilities/Timing.h"
#include "Utilities/FileSystem.h"
#include "GL/freeglut_ext.h"
-#include "Demos/Visualization/Shader.h"
+#include "Visualization/Shader.h"
#include "Simulation/TimeManager.h"
#include "Simulation/Simulation.h"
-#include "Demos/Visualization/Visualization.h"
+#include "Visualization/Visualization.h"
#include "Simulation/TimeStepController.h"
diff --git a/Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDWrapper.h b/Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDWrapper.h
similarity index 100%
rename from Demos/DynamicBoundaryDemo/PositionBasedDynamicsWrapper/PBDWrapper.h
rename to Simulators/DynamicBoundarySimulator/PositionBasedDynamicsWrapper/PBDWrapper.h
diff --git a/Demos/DynamicBoundaryDemo/main.cpp b/Simulators/DynamicBoundarySimulator/main.cpp
similarity index 76%
rename from Demos/DynamicBoundaryDemo/main.cpp
rename to Simulators/DynamicBoundarySimulator/main.cpp
index dc0d3ea..47b99ca 100644
--- a/Demos/DynamicBoundaryDemo/main.cpp
+++ b/Simulators/DynamicBoundarySimulator/main.cpp
@@ -11,13 +11,13 @@
#include "Utilities/OBJLoader.h"
#include "SPlisHSPlasH/Utilities/PoissonDiskSampling.h"
#include "PositionBasedDynamicsWrapper/PBDWrapper.h"
-#include "Demos/Common/DemoBase.h"
+#include "Simulators/Common/SimulatorBase.h"
#include "Utilities/FileSystem.h"
#include "Utilities/Version.h"
#include "Utilities/SystemInfo.h"
#include "Utilities/Counting.h"
#include "SPlisHSPlasH/Simulation.h"
-#include "Demos/Common/TweakBarParameters.h"
+#include "Simulators/Common/TweakBarParameters.h"
#include "SPlisHSPlasH/TimeStep.h"
#include "GL/freeglut_ext.h"
@@ -44,18 +44,27 @@ void updateBoundaryParticles(const bool forceUpdate);
void updateBoundaryForces();
void TW_CALL setCurrentFluidModel(const void *value, void *clientData);
void TW_CALL getCurrentFluidModel(void *value, void *clientData);
-
-DemoBase *base;
+void TW_CALL setColorField(const void *value, void *clientData);
+void TW_CALL getColorField(void *value, void *clientData);
+void TW_CALL setRenderMaxValue(const void *value, void *clientData);
+void TW_CALL getRenderMaxValue(void *value, void *clientData);
+void TW_CALL setRenderMinValue(const void *value, void *clientData);
+void TW_CALL getRenderMinValue(void *value, void *clientData);
+void TW_CALL setColorMapType(const void *value, void *clientData);
+void TW_CALL getColorMapType(void *value, void *clientData);
+
+SimulatorBase *base;
PBDWrapper pbdWrapper;
unsigned int currentFluidModel = 0;
+std::vector<std::string> colorFieldNames;
// main
int main( int argc, char **argv )
{
REPORT_MEMORY_LEAKS;
- base = new DemoBase();
- base->init(argc, argv, "DynamicBoundaryDemo");
+ base = new SimulatorBase();
+ base->init(argc, argv, "DynamicBoundarySimulator");
Simulation *sim = Simulation::getCurrent();
sim->init(base->getScene().particleRadius, base->getScene().sim2D);
@@ -83,6 +92,7 @@ int main( int argc, char **argv )
model->setSurfaceMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getSurfaceTensionBase()); });
model->setViscosityMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getViscosityBase()); });
model->setVorticityMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getVorticityBase()); });
+ model->setElasticityMethodChangedCallback([&]() { reset(); initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getElasticityBase()); });
}
initParameters();
@@ -141,16 +151,46 @@ void initParameters()
// show GUI only for currently selected fluid model
unsigned int i = currentFluidModel;
FluidModel *model = sim->getFluidModel(currentFluidModel);
+
+ colorFieldNames.clear();
+ colorFieldNames.resize(model->numberOfFields());
+ TwType enumType = TwDefineEnum("ColorFieldEnum", NULL, 0);
+ std::ostringstream oss;
+ for (unsigned int j = 0; j < model->numberOfFields(); j++)
+ {
+ const FieldDescription &field = model->getField(j);
+ if ((field.type == FieldType::Scalar) || (field.type == FieldType::Vector3))
+ {
+ if (j != 0)
+ oss << ", ";
+ oss << j << " {" << field.name.c_str() << "}";
+ colorFieldNames[j] = field.name;
+ }
+ }
+ std::string enumStr = " label='Color field' enum='" + oss.str() + "' group='" + model->getId() + "'";
+ enumStr = enumStr + " help='Choose vector or scalar field for particle coloring.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "ColorField", enumType, setColorField, getColorField, nullptr, enumStr.c_str());
+
+ TwType enumType2 = TwDefineEnum("ColorMapTypeEnum", NULL, 0);
+ std::string str = " label='Color map' enum='0 {None}, 1 {Jet}, 2 {Plasma}' group='" + model->getId() + "' help='Choose a color map.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "ColorMapType", enumType2, setColorMapType, getColorMapType, nullptr, str.c_str());
+ str = " label='Min. value (shader)' step=0.001 precision=3 group='" + model->getId() + "' help='Minimal value used for color-coding the color field in the rendering process.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "RenderMinValue", TW_TYPE_REAL, setRenderMinValue, getRenderMinValue, nullptr, str.c_str());
+ str = " label='Max. value (shader)' step=0.001 precision=3 group='" + model->getId() + "' help='Maximal value used for color-coding the color field in the rendering process.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "RenderMaxValue", TW_TYPE_REAL, setRenderMaxValue, getRenderMaxValue, nullptr, str.c_str());
+
TweakBarParameters::createParameterObjectGUI(model);
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getDragBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getSurfaceTensionBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getViscosityBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getVorticityBase());
+ TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getElasticityBase());
TwDefine((std::string("TweakBar/FluidModel group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/'Drag force' group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/'Surface tension' group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/Viscosity group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/Vorticity group='") + model->getId() + "'").c_str());
+ TwDefine((std::string("TweakBar/'Elasticity' group='") + model->getId() + "'").c_str());
pbdWrapper.initGUI();
}
@@ -177,21 +217,21 @@ void reset()
void timeStep ()
{
- const Real stopAt = base->getValue<Real>(DemoBase::STOP_AT);
+ const Real stopAt = base->getValue<Real>(SimulatorBase::STOP_AT);
if ((stopAt > 0.0) && (stopAt < TimeManager::getCurrent()->getTime()))
glutLeaveMainLoop();
- const Real pauseAt = base->getValue<Real>(DemoBase::PAUSE_AT);
+ const Real pauseAt = base->getValue<Real>(SimulatorBase::PAUSE_AT);
if ((pauseAt > 0.0) && (pauseAt < TimeManager::getCurrent()->getTime()))
- base->setValue(DemoBase::PAUSE, true);
+ base->setValue(SimulatorBase::PAUSE, true);
- if (base->getValue<bool>(DemoBase::PAUSE))
+ if (base->getValue<bool>(SimulatorBase::PAUSE))
return;
// Simulation code
Simulation *sim = Simulation::getCurrent();
const bool sim2D = sim->is2DSimulation();
- const unsigned int numSteps = base->getValue<unsigned int>(DemoBase::NUM_STEPS_PER_RENDER);
+ const unsigned int numSteps = base->getValue<unsigned int>(SimulatorBase::NUM_STEPS_PER_RENDER);
for (unsigned int i = 0; i < numSteps; i++)
{
START_TIMING("SimStep");
@@ -235,7 +275,7 @@ void renderBoundary()
Shader &shader = base->getShaderScalar();
Shader &meshShader = base->getMeshShader();
SceneLoader::Scene &scene = base->getScene();
- const int renderWalls = base->getValue<int>(DemoBase::RENDER_WALLS);
+ const int renderWalls = base->getValue<int>(SimulatorBase::RENDER_WALLS);
GLint context_major_version = base->getContextMajorVersion();
if ((renderWalls == 1) || (renderWalls == 2))
@@ -286,7 +326,11 @@ void renderBoundary()
void render()
{
float gridColor[4] = { 0.2f, 0.2f, 0.2f, 1.0f };
- MiniGL::drawGrid(gridColor);
+ const bool sim2D = Simulation::getCurrent()->is2DSimulation();
+ if (sim2D)
+ MiniGL::drawGrid_xy(gridColor);
+ else
+ MiniGL::drawGrid_xz(gridColor);
MiniGL::coordinateSystem();
MiniGL::drawTime(TimeManager::getCurrent()->getTime());
@@ -294,11 +338,9 @@ void render()
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i++)
{
- FluidModel *model = sim->getFluidModel(i);
-
float fluidColor[4] = { 0.3f, 0.5f, 0.9f, 1.0f };
MiniGL::hsvToRgb(0.61f - 0.1f*i, 0.66f, 0.9f, fluidColor);
- base->renderFluid(model, fluidColor);
+ base->renderFluid(i, fluidColor);
}
renderBoundary();
@@ -309,7 +351,7 @@ void render()
PBD::SimulationModel &model = pbdWrapper.getSimulationModel();
PBD::SimulationModel::RigidBodyVector &rb = model.getRigidBodies();
- const int renderWalls = base->getValue<int>(DemoBase::RENDER_WALLS);
+ const int renderWalls = base->getValue<int>(SimulatorBase::RENDER_WALLS);
SceneLoader::Scene &scene = base->getScene();
if ((renderWalls == 3) || (renderWalls == 4))
{
@@ -425,7 +467,7 @@ void initBoundaryData()
}
Simulation::getCurrent()->addBoundaryModel(rb, static_cast<unsigned int>(boundaryParticles.size()), &boundaryParticles[0]);
}
- Simulation::getCurrent()->updateBoundaryPsi();
+ Simulation::getCurrent()->updateBoundaryVolume();
updateBoundaryParticles(true);
}
@@ -488,3 +530,57 @@ void TW_CALL getCurrentFluidModel(void *value, void *clientData)
{
*(unsigned int *)(value) = *((unsigned int*)clientData);
}
+
+void TW_CALL setColorField(const void *value, void *clientData)
+{
+ const unsigned int val = *(const unsigned int *)(value);
+ base->setColorField(currentFluidModel, colorFieldNames[val]);
+}
+
+void TW_CALL getColorField(void *value, void *clientData)
+{
+ const std::string &fieldName = base->getColorField(currentFluidModel);
+ unsigned int index = 0;
+ for (auto i = 0; i < colorFieldNames.size(); i++)
+ {
+ if (colorFieldNames[i] == fieldName)
+ {
+ index = i;
+ break;
+ }
+ }
+ *(unsigned int *)(value) = index;
+}
+
+void TW_CALL setRenderMaxValue(const void *value, void *clientData)
+{
+ const Real val = *(const Real *)(value);
+ base->setRenderMaxValue(currentFluidModel, val);
+}
+
+void TW_CALL getRenderMaxValue(void *value, void *clientData)
+{
+ *(Real *)(value) = base->getRenderMaxValue(currentFluidModel);
+}
+
+void TW_CALL setRenderMinValue(const void *value, void *clientData)
+{
+ const Real val = *(const Real *)(value);
+ base->setRenderMinValue(currentFluidModel, val);
+}
+
+void TW_CALL getRenderMinValue(void *value, void *clientData)
+{
+ *(Real *)(value) = base->getRenderMinValue(currentFluidModel);
+}
+
+void TW_CALL setColorMapType(const void *value, void *clientData)
+{
+ const unsigned int val = *(const unsigned int *)(value);
+ base->setColorMapType(currentFluidModel, val);
+}
+
+void TW_CALL getColorMapType(void *value, void *clientData)
+{
+ *(unsigned int *)(value) = base->getColorMapType(currentFluidModel);
+}
diff --git a/Demos/StaticBoundaryDemo/CMakeLists.txt b/Simulators/StaticBoundarySimulator/CMakeLists.txt
similarity index 66%
rename from Demos/StaticBoundaryDemo/CMakeLists.txt
rename to Simulators/StaticBoundarySimulator/CMakeLists.txt
index 0b00f8c..c1e597e 100644
--- a/Demos/StaticBoundaryDemo/CMakeLists.txt
+++ b/Simulators/StaticBoundarySimulator/CMakeLists.txt
@@ -35,13 +35,13 @@ include_directories(${PROJECT_PATH}/extern/install/GenericParameters/include)
set(SIMULATION_DEPENDENCIES ${SIMULATION_DEPENDENCIES} Ext_GenericParameters)
-add_executable(StaticBoundaryDemo
+add_executable(StaticBoundarySimulator
main.cpp
- ${PROJECT_PATH}/Demos/Common/DemoBase.cpp
- ${PROJECT_PATH}/Demos/Common/DemoBase.h
- ${PROJECT_PATH}/Demos/Common/TweakBarParameters.cpp
- ${PROJECT_PATH}/Demos/Common/TweakBarParameters.h
+ ${PROJECT_PATH}/Simulators/Common/SimulatorBase.cpp
+ ${PROJECT_PATH}/Simulators/Common/SimulatorBase.h
+ ${PROJECT_PATH}/Simulators/Common/TweakBarParameters.cpp
+ ${PROJECT_PATH}/Simulators/Common/TweakBarParameters.h
${VIS_FILES}
)
@@ -55,15 +55,15 @@ add_definitions(-DTW_NO_LIB_PRAGMA -DTW_STATIC)
include_directories(${PROJECT_PATH}/extern/freeglut/include)
include_directories(${PROJECT_PATH}/extern/glew/include)
-set_target_properties(StaticBoundaryDemo PROPERTIES DEBUG_POSTFIX ${CMAKE_DEBUG_POSTFIX})
-set_target_properties(StaticBoundaryDemo PROPERTIES RELWITHDEBINFO_POSTFIX ${CMAKE_RELWITHDEBINFO_POSTFIX})
-set_target_properties(StaticBoundaryDemo PROPERTIES MINSIZEREL_POSTFIX ${CMAKE_MINSIZEREL_POSTFIX})
-add_dependencies(StaticBoundaryDemo ${SIMULATION_DEPENDENCIES})
-target_link_libraries(StaticBoundaryDemo ${SIMULATION_LINK_LIBRARIES})
+set_target_properties(StaticBoundarySimulator PROPERTIES DEBUG_POSTFIX ${CMAKE_DEBUG_POSTFIX})
+set_target_properties(StaticBoundarySimulator PROPERTIES RELWITHDEBINFO_POSTFIX ${CMAKE_RELWITHDEBINFO_POSTFIX})
+set_target_properties(StaticBoundarySimulator PROPERTIES MINSIZEREL_POSTFIX ${CMAKE_MINSIZEREL_POSTFIX})
+add_dependencies(StaticBoundarySimulator ${SIMULATION_DEPENDENCIES})
+target_link_libraries(StaticBoundarySimulator ${SIMULATION_LINK_LIBRARIES})
VIS_SOURCE_GROUPS()
-set_target_properties(StaticBoundaryDemo PROPERTIES FOLDER "Demos")
+set_target_properties(StaticBoundarySimulator PROPERTIES FOLDER "Simulators")
-add_custom_command(TARGET StaticBoundaryDemo PRE_BUILD
+add_custom_command(TARGET StaticBoundarySimulator PRE_BUILD
COMMAND ${CMAKE_COMMAND} -E copy_directory
- ${PROJECT_PATH}/data/shaders $<TARGET_FILE_DIR:StaticBoundaryDemo>/resources/shaders)
+ ${PROJECT_PATH}/data/shaders $<TARGET_FILE_DIR:StaticBoundarySimulator>/resources/shaders)
diff --git a/Demos/StaticBoundaryDemo/main.cpp b/Simulators/StaticBoundarySimulator/main.cpp
similarity index 74%
rename from Demos/StaticBoundaryDemo/main.cpp
rename to Simulators/StaticBoundarySimulator/main.cpp
index eddcb6e..1670585 100644
--- a/Demos/StaticBoundaryDemo/main.cpp
+++ b/Simulators/StaticBoundarySimulator/main.cpp
@@ -10,14 +10,14 @@
#include "SPlisHSPlasH/StaticRigidBody.h"
#include "Utilities/OBJLoader.h"
#include "SPlisHSPlasH/Utilities/PoissonDiskSampling.h"
-#include "Demos/Common/DemoBase.h"
+#include "Simulators/Common/SimulatorBase.h"
#include "Utilities/FileSystem.h"
#include "Utilities/Version.h"
#include <fstream>
#include "Utilities/Logger.h"
#include "Utilities/Counting.h"
#include "SPlisHSPlasH/Simulation.h"
-#include "Demos/Common/TweakBarParameters.h"
+#include "Simulators/Common/TweakBarParameters.h"
#include "GL/freeglut_ext.h"
// Enable memory leak detection
@@ -42,17 +42,26 @@ void initParameters();
void loadObj(const std::string &filename, TriangleMesh &geo, const Vector3r &scale);
void TW_CALL setCurrentFluidModel(const void *value, void *clientData);
void TW_CALL getCurrentFluidModel(void *value, void *clientData);
-
-DemoBase *base;
+void TW_CALL setColorField(const void *value, void *clientData);
+void TW_CALL getColorField(void *value, void *clientData);
+void TW_CALL setRenderMaxValue(const void *value, void *clientData);
+void TW_CALL getRenderMaxValue(void *value, void *clientData);
+void TW_CALL setRenderMinValue(const void *value, void *clientData);
+void TW_CALL getRenderMinValue(void *value, void *clientData);
+void TW_CALL setColorMapType(const void *value, void *clientData);
+void TW_CALL getColorMapType(void *value, void *clientData);
+
+SimulatorBase *base;
unsigned int currentFluidModel = 0;
+std::vector<std::string> colorFieldNames;
// main
int main( int argc, char **argv )
{
REPORT_MEMORY_LEAKS;
- base = new DemoBase();
- base->init(argc, argv, "StaticBoundaryDemo");
+ base = new SimulatorBase();
+ base->init(argc, argv, "StaticBoundarySimulator");
Simulation *sim = Simulation::getCurrent();
sim->init(base->getScene().particleRadius, base->getScene().sim2D);
@@ -74,6 +83,7 @@ int main( int argc, char **argv )
model->setSurfaceMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getSurfaceTensionBase()); });
model->setViscosityMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getViscosityBase()); });
model->setVorticityMethodChangedCallback([&]() { initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getVorticityBase()); });
+ model->setElasticityMethodChangedCallback([&]() { reset(); initParameters(); base->getSceneLoader()->readParameterObject(key, (ParameterObject*)model->getElasticityBase()); });
}
initParameters();
@@ -130,16 +140,46 @@ void initParameters()
// show GUI only for currently selected fluid model
unsigned int i = currentFluidModel;
FluidModel *model = sim->getFluidModel(currentFluidModel);
+
+ colorFieldNames.clear();
+ colorFieldNames.resize(model->numberOfFields());
+ TwType enumType = TwDefineEnum("ColorFieldEnum", NULL, 0);
+ std::ostringstream oss;
+ for (unsigned int j = 0; j < model->numberOfFields(); j++)
+ {
+ const FieldDescription &field = model->getField(j);
+ if ((field.type == FieldType::Scalar) || (field.type == FieldType::Vector3))
+ {
+ if (j != 0)
+ oss << ", ";
+ oss << j << " {" << field.name.c_str() << "}";
+ colorFieldNames[j] = field.name;
+ }
+ }
+ std::string enumStr = " label='Color field' enum='" + oss.str() + "' group='" + model->getId() + "'";
+ enumStr = enumStr + " help='Choose vector or scalar field for particle coloring.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "ColorField", enumType, setColorField, getColorField, nullptr, enumStr.c_str());
+
+ TwType enumType2 = TwDefineEnum("ColorMapTypeEnum", NULL, 0);
+ std::string str = " label='Color map' enum='0 {None}, 1 {Jet}, 2 {Plasma}' group='" + model->getId() + "' help='Choose a color map.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "ColorMapType", enumType2, setColorMapType, getColorMapType, nullptr, str.c_str());
+ str = " label='Min. value (shader)' step=0.001 precision=3 group='" + model->getId() + "' help='Minimal value used for color-coding the color field in the rendering process.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "RenderMinValue", TW_TYPE_REAL, setRenderMinValue, getRenderMinValue, nullptr, str.c_str());
+ str = " label='Max. value (shader)' step=0.001 precision=3 group='" + model->getId() + "' help='Maximal value used for color-coding the color field in the rendering process.'";
+ TwAddVarCB(MiniGL::getTweakBar(), "RenderMaxValue", TW_TYPE_REAL, setRenderMaxValue, getRenderMaxValue, nullptr, str.c_str());
+
TweakBarParameters::createParameterObjectGUI(model);
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getDragBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getSurfaceTensionBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getViscosityBase());
TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getVorticityBase());
+ TweakBarParameters::createParameterObjectGUI((GenParam::ParameterObject*) model->getElasticityBase());
TwDefine((std::string("TweakBar/FluidModel group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/'Drag force' group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/'Surface tension' group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/Viscosity group='") + model->getId() + "'").c_str());
TwDefine((std::string("TweakBar/Vorticity group='") + model->getId() + "'").c_str());
+ TwDefine((std::string("TweakBar/'Elasticity' group='") + model->getId() + "'").c_str());
}
void reset()
@@ -157,21 +197,21 @@ void reset()
void timeStep ()
{
- const Real stopAt = base->getValue<Real>(DemoBase::STOP_AT);
+ const Real stopAt = base->getValue<Real>(SimulatorBase::STOP_AT);
if ((stopAt > 0.0) && (stopAt < TimeManager::getCurrent()->getTime()))
glutLeaveMainLoop();
- const Real pauseAt = base->getValue<Real>(DemoBase::PAUSE_AT);
+ const Real pauseAt = base->getValue<Real>(SimulatorBase::PAUSE_AT);
if ((pauseAt > 0.0) && (pauseAt < TimeManager::getCurrent()->getTime()))
- base->setValue(DemoBase::PAUSE, true);
+ base->setValue(SimulatorBase::PAUSE, true);
- if (base->getValue<bool>(DemoBase::PAUSE))
+ if (base->getValue<bool>(SimulatorBase::PAUSE))
return;
// Simulation code
Simulation *sim = Simulation::getCurrent();
const bool sim2D = sim->is2DSimulation();
- const unsigned int numSteps = base->getValue<unsigned int>(DemoBase::NUM_STEPS_PER_RENDER);
+ const unsigned int numSteps = base->getValue<unsigned int>(SimulatorBase::NUM_STEPS_PER_RENDER);
for (unsigned int i = 0; i < numSteps; i++)
{
START_TIMING("SimStep");
@@ -201,7 +241,11 @@ void timeStep ()
void render()
{
float gridColor[4] = { 0.2f, 0.2f, 0.2f, 1.0f };
- MiniGL::drawGrid(gridColor);
+ const bool sim2D = Simulation::getCurrent()->is2DSimulation();
+ if (sim2D)
+ MiniGL::drawGrid_xy(gridColor);
+ else
+ MiniGL::drawGrid_xz(gridColor);
MiniGL::coordinateSystem();
MiniGL::drawTime(TimeManager::getCurrent()->getTime());
@@ -209,11 +253,9 @@ void render()
Simulation *sim = Simulation::getCurrent();
for (unsigned int i = 0; i < sim->numberOfFluidModels(); i++)
{
- FluidModel *model = sim->getFluidModel(i);
-
float fluidColor[4] = { 0.3f, 0.5f, 0.9f, 1.0f };
MiniGL::hsvToRgb(0.61f-0.1f*i, 0.66f, 0.9f, fluidColor);
- base->renderFluid(model, fluidColor);
+ base->renderFluid(i, fluidColor);
}
renderBoundary();
}
@@ -224,7 +266,7 @@ void renderBoundary()
Shader &shader = base->getShaderScalar();
Shader &meshShader = base->getMeshShader();
SceneLoader::Scene &scene = base->getScene();
- const int renderWalls = base->getValue<int>(DemoBase::RENDER_WALLS);
+ const int renderWalls = base->getValue<int>(SimulatorBase::RENDER_WALLS);
GLint context_major_version = base->getContextMajorVersion();
if ((renderWalls == 1) || (renderWalls == 2))
@@ -412,7 +454,7 @@ void initBoundaryData()
Simulation::getCurrent()->addBoundaryModel(rb, static_cast<unsigned int>(boundaryParticles.size()), &boundaryParticles[0]);
}
- Simulation::getCurrent()->updateBoundaryPsi();
+ Simulation::getCurrent()->updateBoundaryVolume();
}
void TW_CALL setCurrentFluidModel(const void *value, void *clientData)
@@ -427,3 +469,57 @@ void TW_CALL getCurrentFluidModel(void *value, void *clientData)
{
*(unsigned int *)(value) = *((unsigned int*)clientData);
}
+
+void TW_CALL setColorField(const void *value, void *clientData)
+{
+ const unsigned int val = *(const unsigned int *)(value);
+ base->setColorField(currentFluidModel, colorFieldNames[val]);
+}
+
+void TW_CALL getColorField(void *value, void *clientData)
+{
+ const std::string &fieldName = base->getColorField(currentFluidModel);
+ unsigned int index = 0;
+ for (auto i = 0; i < colorFieldNames.size(); i++)
+ {
+ if (colorFieldNames[i] == fieldName)
+ {
+ index = i;
+ break;
+ }
+ }
+ *(unsigned int *)(value) = index;
+}
+
+void TW_CALL setRenderMaxValue(const void *value, void *clientData)
+{
+ const Real val = *(const Real *)(value);
+ base->setRenderMaxValue(currentFluidModel, val);
+}
+
+void TW_CALL getRenderMaxValue(void *value, void *clientData)
+{
+ *(Real *)(value) = base->getRenderMaxValue(currentFluidModel);
+}
+
+void TW_CALL setRenderMinValue(const void *value, void *clientData)
+{
+ const Real val = *(const Real *)(value);
+ base->setRenderMinValue(currentFluidModel, val);
+}
+
+void TW_CALL getRenderMinValue(void *value, void *clientData)
+{
+ *(Real *)(value) = base->getRenderMinValue(currentFluidModel);
+}
+
+void TW_CALL setColorMapType(const void *value, void *clientData)
+{
+ const unsigned int val = *(const unsigned int *)(value);
+ base->setColorMapType(currentFluidModel, val);
+}
+
+void TW_CALL getColorMapType(void *value, void *clientData)
+{
+ *(unsigned int *)(value) = base->getColorMapType(currentFluidModel);
+}
diff --git a/Visualization/MiniGL.cpp b/Visualization/MiniGL.cpp
index 31b7147..f9580e6 100644
--- a/Visualization/MiniGL.cpp
+++ b/Visualization/MiniGL.cpp
@@ -393,7 +393,7 @@ void MiniGL::drawTetrahedron(const Vector3r &a, const Vector3r &b, const Vector3
drawTriangle(b, c, d, normal4, color);
}
-void MiniGL::drawGrid(float *color)
+void MiniGL::drawGrid_xz(float *color)
{
float speccolor[4] = { 1.0, 1.0, 1.0, 1.0 };
glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, color);
@@ -422,6 +422,35 @@ void MiniGL::drawGrid(float *color)
glEnd();
}
+void MiniGL::drawGrid_xy(float *color)
+{
+ float speccolor[4] = { 1.0, 1.0, 1.0, 1.0 };
+ glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT, color);
+ glMaterialfv(GL_FRONT_AND_BACK, GL_DIFFUSE, color);
+ glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, speccolor);
+ glMaterialf(GL_FRONT_AND_BACK, GL_SHININESS, 100.0);
+
+ const int size = 5;
+
+ glBegin(GL_LINES);
+ for (int i = -size; i <= size; i++)
+ {
+ glVertex3f((float)i, (float)-size, 0.0f );
+ glVertex3f((float)i, (float)size, 0.0f);
+ glVertex3f((float)-size, (float)i, 0.0f);
+ glVertex3f((float)size, (float)i, 0.0f);
+ }
+ glEnd();
+
+ glLineWidth(3.0f);
+ glBegin(GL_LINES);
+ glVertex3f((float)-size, 0.0f, 0.0f);
+ glVertex3f((float)size, 0.0f, 0.0f);
+ glVertex3f(0.0f, (float)-size, 0.0f);
+ glVertex3f(0.0f, (float)size, 0.0f);
+ glEnd();
+}
+
void MiniGL::setViewport(float pfovy, float pznear, float pzfar, const Vector3r &peyepoint, const Vector3r &plookat)
{
fovy = pfovy;
diff --git a/Visualization/MiniGL.h b/Visualization/MiniGL.h
index 70bcdde..2f65fd2 100644
--- a/Visualization/MiniGL.h
+++ b/Visualization/MiniGL.h
@@ -129,7 +129,8 @@ namespace SPH
*/
static void drawTetrahedron(const Vector3r &a, const Vector3r &b, const Vector3r &c, const Vector3r &d, float *color);
static void drawTriangle (const Vector3r &a, const Vector3r &b, const Vector3r &c, const Vector3r &norm, float *color);
- static void drawGrid(float *color);
+ static void drawGrid_xz(float *color);
+ static void drawGrid_xy(float *color);
static void drawBitmapText (float x, float y, const char *str, int strLength, float *color);
static void drawStrokeText(const Real x, const Real y, const Real z, float scale, const char *str, int strLength, float *color);
static void drawStrokeText (const Vector3r &pos, float scale, const char *str, int strLength, float *color);
diff --git a/bin/Buckling.bat b/bin/Buckling.bat
deleted file mode 100644
index 5b028a7..0000000
--- a/bin/Buckling.bat
+++ /dev/null
@@ -1,2 +0,0 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel.json
-pause
\ No newline at end of file
diff --git a/bin/Buckling_Bender2017.bat b/bin/Buckling_Bender2017.bat
index 61f762a..44b5c2d 100644
--- a/bin/Buckling_Bender2017.bat
+++ b/bin/Buckling_Bender2017.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel_Bender2017.json
+.\StaticBoundarySimulator.exe ../data/scenes/BucklingModel_Bender2017.json
pause
\ No newline at end of file
diff --git a/bin/Buckling_Peer2015.bat b/bin/Buckling_Peer2015.bat
index 95a9e75..1fedbfc 100644
--- a/bin/Buckling_Peer2015.bat
+++ b/bin/Buckling_Peer2015.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel_Peer2015.json
+.\StaticBoundarySimulator.exe ../data/scenes/BucklingModel_Peer2015.json
pause
\ No newline at end of file
diff --git a/bin/Buckling_Peer2016.bat b/bin/Buckling_Peer2016.bat
index 4734a82..1ede812 100644
--- a/bin/Buckling_Peer2016.bat
+++ b/bin/Buckling_Peer2016.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel_Peer2016.json
+.\StaticBoundarySimulator.exe ../data/scenes/BucklingModel_Peer2016.json
pause
\ No newline at end of file
diff --git a/bin/Buckling_Takahashi2015.bat b/bin/Buckling_Takahashi2015.bat
index e6cb2c4..f439ae6 100644
--- a/bin/Buckling_Takahashi2015.bat
+++ b/bin/Buckling_Takahashi2015.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel_Takahashi2015.json
+.\StaticBoundarySimulator.exe ../data/scenes/BucklingModel_Takahashi2015.json
pause
\ No newline at end of file
diff --git a/bin/Buckling_Weiler2018.bat b/bin/Buckling_Weiler2018.bat
index 26b5d1d..d314731 100644
--- a/bin/Buckling_Weiler2018.bat
+++ b/bin/Buckling_Weiler2018.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/BucklingModel_Weiler2018.json
+.\StaticBoundarySimulator.exe ../data/scenes/BucklingModel_Weiler2018.json
pause
\ No newline at end of file
diff --git a/bin/Bunny_vs_Dragon.bat b/bin/Bunny_vs_Dragon.bat
index c995168..cfce186 100644
--- a/bin/Bunny_vs_Dragon.bat
+++ b/bin/Bunny_vs_Dragon.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/Bunny_vs_Dragon.json
+.\StaticBoundarySimulator.exe ../data/scenes/Bunny_vs_Dragon.json
pause
\ No newline at end of file
diff --git a/bin/Coiling_Bender2017.bat b/bin/Coiling_Bender2017.bat
index 489785c..cf2cc7a 100644
--- a/bin/Coiling_Bender2017.bat
+++ b/bin/Coiling_Bender2017.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/CoilingModel_Bender2017.json
+.\StaticBoundarySimulator.exe ../data/scenes/CoilingModel_Bender2017.json
pause
\ No newline at end of file
diff --git a/bin/Coiling_Peer2015.bat b/bin/Coiling_Peer2015.bat
index 467569c..a78a454 100644
--- a/bin/Coiling_Peer2015.bat
+++ b/bin/Coiling_Peer2015.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/CoilingModel_Peer2015.json
+.\StaticBoundarySimulator.exe ../data/scenes/CoilingModel_Peer2015.json
pause
\ No newline at end of file
diff --git a/bin/Coiling_Peer2016.bat b/bin/Coiling_Peer2016.bat
index b851695..2523f2a 100644
--- a/bin/Coiling_Peer2016.bat
+++ b/bin/Coiling_Peer2016.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/CoilingModel_Peer2016.json
+.\StaticBoundarySimulator.exe ../data/scenes/CoilingModel_Peer2016.json
pause
\ No newline at end of file
diff --git a/bin/Coiling_Takahashi2015.bat b/bin/Coiling_Takahashi2015.bat
index 3c44b1d..23f88f2 100644
--- a/bin/Coiling_Takahashi2015.bat
+++ b/bin/Coiling_Takahashi2015.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/CoilingModel_Takahashi2015.json
+.\StaticBoundarySimulator.exe ../data/scenes/CoilingModel_Takahashi2015.json
pause
\ No newline at end of file
diff --git a/bin/Coiling_Weiler2018.bat b/bin/Coiling_Weiler2018.bat
index c5967c5..0242a0d 100644
--- a/bin/Coiling_Weiler2018.bat
+++ b/bin/Coiling_Weiler2018.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/CoilingModel_Weiler2018.json
+.\StaticBoundarySimulator.exe ../data/scenes/CoilingModel_Weiler2018.json
pause
\ No newline at end of file
diff --git a/bin/DamBreakModel.bat b/bin/DamBreakModel.bat
index 488e258..d0ad0b2 100644
--- a/bin/DamBreakModel.bat
+++ b/bin/DamBreakModel.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DamBreakModel.json
+.\StaticBoundarySimulator.exe ../data/scenes/DamBreakModel.json
pause
\ No newline at end of file
diff --git a/bin/DamBreakModelDragons.bat b/bin/DamBreakModelDragons.bat
index c2e81dd..feba9c0 100644
--- a/bin/DamBreakModelDragons.bat
+++ b/bin/DamBreakModelDragons.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DamBreakModelDragons.json
+.\StaticBoundarySimulator.exe ../data/scenes/DamBreakModelDragons.json
pause
\ No newline at end of file
diff --git a/bin/DamBreakModel_2D.bat b/bin/DamBreakModel_2D.bat
index f110733..bc9163a 100644
--- a/bin/DamBreakModel_2D.bat
+++ b/bin/DamBreakModel_2D.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DamBreakModel_2D.json
+.\StaticBoundarySimulator.exe ../data/scenes/DamBreakModel_2D.json
pause
\ No newline at end of file
diff --git a/bin/DamBreakWithObjects.bat b/bin/DamBreakWithObjects.bat
index 91913cf..853eb2b 100644
--- a/bin/DamBreakWithObjects.bat
+++ b/bin/DamBreakWithObjects.bat
@@ -1,2 +1,2 @@
-.\DynamicBoundaryDemo.exe ../data/scenes/DamBreakWithObjects.json
+.\DynamicBoundarySimulator.exe ../data/scenes/DamBreakWithObjects.json
pause
\ No newline at end of file
diff --git a/bin/DoubleDamBreak.bat b/bin/DoubleDamBreak.bat
index f834a05..c1adff8 100644
--- a/bin/DoubleDamBreak.bat
+++ b/bin/DoubleDamBreak.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DoubleDamBreak.json
+.\StaticBoundarySimulator.exe ../data/scenes/DoubleDamBreak.json
pause
\ No newline at end of file
diff --git a/bin/DoubleDamBreakMultiPhase.bat b/bin/DoubleDamBreakMultiPhase.bat
new file mode 100644
index 0000000..16c8658
--- /dev/null
+++ b/bin/DoubleDamBreakMultiPhase.bat
@@ -0,0 +1,2 @@
+.\StaticBoundarySimulator.exe ../data/scenes/DoubleDamBreakMultiPhase.json
+pause
\ No newline at end of file
diff --git a/bin/DoubleDamBreakTwoPhase.bat b/bin/DoubleDamBreakTwoPhase.bat
deleted file mode 100644
index a8fb926..0000000
--- a/bin/DoubleDamBreakTwoPhase.bat
+++ /dev/null
@@ -1,2 +0,0 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DoubleDamBreakTwoPhase.json
-pause
\ No newline at end of file
diff --git a/bin/DoubleDamBreakWithSphere.bat b/bin/DoubleDamBreakWithSphere.bat
index b9d5833..0fec86e 100644
--- a/bin/DoubleDamBreakWithSphere.bat
+++ b/bin/DoubleDamBreakWithSphere.bat
@@ -1,2 +1,2 @@
-.\DynamicBoundaryDemo.exe ../data/scenes/DoubleDamBreakWithSphere.json
+.\DynamicBoundarySimulator.exe ../data/scenes/DoubleDamBreakWithSphere.json
pause
\ No newline at end of file
diff --git a/bin/DragTest.bat b/bin/DragTest.bat
index 91b795b..b22289e 100644
--- a/bin/DragTest.bat
+++ b/bin/DragTest.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/DragTest.json
+.\StaticBoundarySimulator.exe ../data/scenes/DragTest.json
pause
diff --git a/bin/Emitter.bat b/bin/Emitter.bat
index 472a7d2..2360991 100644
--- a/bin/Emitter.bat
+++ b/bin/Emitter.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/Emitter.json
+.\StaticBoundarySimulator.exe ../data/scenes/Emitter.json
pause
\ No newline at end of file
diff --git a/bin/MultiPhaseColoring.bat b/bin/MultiPhaseColoring.bat
new file mode 100644
index 0000000..ce96122
--- /dev/null
+++ b/bin/MultiPhaseColoring.bat
@@ -0,0 +1,2 @@
+.\StaticBoundarySimulator.exe ../data/scenes/MultiPhaseColoring.json
+pause
\ No newline at end of file
diff --git a/bin/Obstacle.bat b/bin/Obstacle.bat
index a74854c..5a2e33e 100644
--- a/bin/Obstacle.bat
+++ b/bin/Obstacle.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/Obstacle.json
+.\StaticBoundarySimulator.exe ../data/scenes/Obstacle.json
pause
\ No newline at end of file
diff --git a/bin/ViscousBunny.bat b/bin/ViscousBunny.bat
index 2fcd3d9..0ad070e 100644
--- a/bin/ViscousBunny.bat
+++ b/bin/ViscousBunny.bat
@@ -1,2 +1,2 @@
-.\StaticBoundaryDemo.exe ../data/scenes/ViscousBunny.json
+.\StaticBoundarySimulator.exe ../data/scenes/ViscousBunny.json
pause
\ No newline at end of file
diff --git a/data/Scenes/DamBreakModel.json b/data/Scenes/DamBreakModel.json
index 1236a39..6899760 100644
--- a/data/Scenes/DamBreakModel.json
+++ b/data/Scenes/DamBreakModel.json
@@ -16,7 +16,8 @@
"stiffness": 50000,
"exponent": 7,
"velocityUpdateMethod": 0,
- "enableDivergenceSolver": true
+ "enableDivergenceSolver": true,
+ "partioAttributes": "density;velocity"
},
"Fluid":
{
diff --git a/data/Scenes/DeformableModel.json b/data/Scenes/DeformableModel.json
new file mode 100644
index 0000000..249dbce
--- /dev/null
+++ b/data/Scenes/DeformableModel.json
@@ -0,0 +1,55 @@
+{
+ "Configuration":
+ {
+ "imteStepSize": 0.001,
+ "particleRadius": 0.025,
+ "numberOfStepsPerRenderUpdate": 2,
+ "density0": 1000,
+ "simulationMethod": 4,
+ "gravitation": [0,-9.81,0],
+ "cflMethod": 0,
+ "cflFactor": 1,
+ "cflMaxTimeStepSize": 0.005,
+ "maxIterations": 100,
+ "maxError": 0.01,
+ "maxIterationsV": 100,
+ "maxErrorV": 0.1,
+ "stiffness": 50000,
+ "exponent": 7,
+ "velocityUpdateMethod": 0,
+ "enableDivergenceSolver": true
+ },
+ "Fluid":
+ {
+ "elasticityMethod": 2,
+ "viscosity": 0.01,
+ "viscosityMethod": 1,
+ "surfaceTension": 0.02,
+ "surfaceTensionMethod": 0,
+ "youngsModulus": 250000.0,
+ "poissonsRatio": 0.33,
+ "alpha": 0.1
+ },
+ "RigidBodies": [
+ {
+ "geometryFile": "../models/UnitBox.obj",
+ "translation": [0,2.5,0],
+ "rotationAxis": [1, 0, 0],
+ "rotationAngle": 0,
+ "scale": [3.1, 5, 3.1],
+ "color": [0.1, 0.4, 0.6, 1.0],
+ "isDynamic": false,
+ "isWall": true
+ }
+ ],
+ "FluidBlocks": [
+ {
+ "denseMode": 0,
+ "start": [-0.25, 0.5, -0.25],
+ "end": [0.25, 1.0, 0.25]
+ }
+ ]
+}
+
+
+
diff --git a/data/Scenes/DoubleDamBreakTwoPhase.json b/data/Scenes/DoubleDamBreakMultiPhase.json
similarity index 100%
rename from data/Scenes/DoubleDamBreakTwoPhase.json
rename to data/Scenes/DoubleDamBreakMultiPhase.json
diff --git a/data/Scenes/MultiPhaseColoring.json b/data/Scenes/MultiPhaseColoring.json
new file mode 100644
index 0000000..f70743d
--- /dev/null
+++ b/data/Scenes/MultiPhaseColoring.json
@@ -0,0 +1,124 @@
+{
+ "Configuration":
+ {
+ "timeStepSize": 0.001,
+ "numberOfStepsPerRenderUpdate": 2,
+ "particleRadius": 0.025,
+ "simulationMethod": 4,
+ "gravitation": [0.0,-9.81,0],
+ "cflMethod": 1,
+ "cflFactor": 0.5,
+ "cflMaxTimeStepSize": 0.005,
+ "maxIterations": 100,
+ "maxError": 0.01,
+ "maxIterationsV": 100,
+ "maxErrorV": 0.1,
+ "stiffness": 5000,
+ "exponent": 1,
+ "velocityUpdateMethod": 0,
+ "enableDivergenceSolver": true
+ },
+ "first_fluid":
+ {
+ "density0": 1000,
+ "viscosity": 0.01,
+ "viscosityBoundary": 5000,
+ "viscosityMethod": 1,
+ "viscoMaxIter": 200,
+ "viscoMaxError": 0.05,
+ "maxEmitterParticles": 6000,
+ "colorField": "velocity",
+ "colorMapType": 1,
+ "renderMinValue": 0.01,
+ "renderMaxValue": 6.0
+ },
+ "second_fluid":
+ {
+ "density0": 100,
+ "viscosity": 0.01,
+ "viscosityMethod": 1,
+ "vorticityMethod": 1,
+ "vorticity": 0.1,
+ "colorField": "angular velocity",
+ "colorMapType": 2,
+ "renderMinValue": 0.01,
+ "renderMaxValue": 2.0
+ },
+ "third_fluid":
+ {
+ "density0": 50,
+ "viscosity": 0.015,
+ "maxEmitterParticles": 1000,
+ "emitterReuseParticles": false,
+ "emitterBoxMin": [-4.0,-1.0,-4.0],
+ "emitterBoxMax": [0.0,4,4.0],
+ "colorField": "position",
+ "colorMapType": 1,
+ "renderMinValue": 0.0,
+ "renderMaxValue": 2.0
+ },
+ "RigidBodies": [
+ {
+ "geometryFile": "../models/UnitBox.obj",
+ "translation": [0,2,0],
+ "rotationAxis": [1, 0, 0],
+ "rotationAngle": 0,
+ "scale": [2.1, 4, 2.1],
+ "color": [0.1, 0.4, 0.6, 1.0],
+ "isDynamic": false,
+ "isWall": true,
+ "restitution" : 0.6,
+ "friction" : 0.0,
+ "collisionObjectType": 2,
+ "collisionObjectScale": [2.1, 4, 2.1],
+ "invertSDF": true
+ }
+ ],
+ "Comment": [
+ {
+ "id": 1,
+ "geometryFile": "../models/torus.obj",
+ "isDynamic": 1,
+ "density": 50,
+ "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": 4,
+ "collisionObjectScale": [0.2, 0.1, 0.2]
+ }
+ ],
+ "FluidBlocks": [
+ {
+ "id": "first_fluid",
+ "denseMode": 0,
+ "start": [-1.0, 0.0, -1.0],
+ "end": [-0.3, 0.75, -0.3]
+ },
+ {
+ "id": "second_fluid",
+ "denseMode": 0,
+ "start": [0.3, 0.0, 0.3],
+ "end": [1.0, 1.75, 1.0]
+ }
+ ],
+ "Emitters": [
+ {
+ "id": "third_fluid",
+ "width": 8,
+ "height": 3,
+ "translation": [-0.2,0.75,0.0],
+ "direction": [1, 0, 0],
+ "velocity": [2, 0, 0.0],
+ "comment" : "The initial velocity should be at least 2*particleRadius*emitsPerSecond.",
+ "emitsPerSecond": 40,
+ "type": 0
+ }
+ ]
+}
+
+
+
diff --git a/data/Scenes/Obstacle.json b/data/Scenes/Obstacle.json
index a4e4a9c..3577324 100644
--- a/data/Scenes/Obstacle.json
+++ b/data/Scenes/Obstacle.json
@@ -25,7 +25,7 @@
"vorticityMethod": 1,
"viscosity": 0.01,
"vorticity": 0.15,
- "viscosityOmega": 0.1,
+ "viscosityOmega": 0.5,
"inertiaInverse": 0.5,
"viscosityMethod": 1,
"surfaceTension": 0.02,
--
GitLab