- 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

@@ -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		
 ) 

Changelog.txt

+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

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.

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;
 }
 

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;
-			}
 	};
 }
 

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})
 

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();

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;

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();
 

SPlisHSPlasH/Elasticity/ElasticityBase.cpp

+#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);
+}
+
+

SPlisHSPlasH/Elasticity/ElasticityBase.h

+#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

SPlisHSPlasH/Elasticity/Elasticity_Becker2009.cpp

+#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)
+			{
+				//////////////////////////////////////////////////////////////////////////
+				// Ganzenmller, G.C. 2015. An hourglass control algorithm for Lagrangian 
+				// Smooth Particle Hydrodynamics. Computer Methods in Applied Mechanics and 
+				// Engineering 286, 87106.
+				//////////////////////////////////////////////////////////////////////////
+				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: Ganzenmller defines xij = xj-xi
+					const Vector3r xi_xj = -(xi - xj);
+					const Real xixj_l = xi_xj.norm();
+					if (xixj_l > 1.0e-6)
+					{
+						// Note: Ganzenmller 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);
+		}
+	}
+}
+

SPlisHSPlasH/Elasticity/Elasticity_Becker2009.h

+#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

SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h

+#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

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;
+		}
+	}
+}
+

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,6 +70,9 @@ 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)

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;

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);
 		}
 	}
 }

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)
 	{