- added XPBD FEM constraint for deformable solids

- updated to Eigen 3.4.0

Changelog.txt

@@ -1,3 +1,7 @@
+2.1.0
+	- added XPBD FEM constraint for deformable solids
+	- updated to Eigen 3.4.0
+
 2.0.1
 	- extended and modified Python interface
 	- more Python examples

Demos/BarDemo/main.cpp

@@ -79,7 +79,8 @@ int main( int argc, char **argv )
 
 	TwType enumType2 = TwDefineEnum("SimulationMethodType", NULL, 0);
 	TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType2, setSimulationMethod, getSimulationMethod, &simulationMethod,
-			" label='Simulation method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics (no inversion handling)}, 4 {Shape matching (no inversion handling)}, 5 {XPBD volume constraints}' group=Simulation");
+			" label='Simulation method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {FEM based XPBD}, \
+			4 {Strain based dynamics (no inversion handling)}, 5 {Shape matching (no inversion handling)}, 6 {XPBD volume constraints}' group=Simulation");
 	TwAddVarCB(MiniGL::getTweakBar(), "stiffness", TW_TYPE_REAL, setStiffness, getStiffness, model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Solid' ");
 	TwAddVarCB(MiniGL::getTweakBar(), "poissonRatio", TW_TYPE_REAL, setPoissonRatio, getPoissonRatio, model, " label='Poisson ratio'  min=0.0 step=0.1 precision=4 group='Solid' ");
 	TwAddVarCB(MiniGL::getTweakBar(), "normalizeStretch", TW_TYPE_BOOL32, setNormalizeStretch, getNormalizeStretch, model, " label='Normalize stretch' group='Solid' ");
@@ -167,11 +168,13 @@ void createMesh()
 
 	// init constraints
 	stiffness = 1.0;
-	if (simulationMethod == 5)
+	if (simulationMethod == 3) 
+		stiffness = 1000000;
+	if (simulationMethod == 6)
 		stiffness = 100000;
 
 	volumeStiffness = 1.0;
-	if (simulationMethod == 5)
+	if (simulationMethod == 6)
 		volumeStiffness = 100000;
 	for (unsigned int cm = 0; cm < model->getTetModels().size(); cm++)
 	{
@@ -201,6 +204,7 @@ void TW_CALL setStiffness(const void* value, void* clientData)
 {
 	stiffness = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_stiffness>(stiffness);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_stiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_stretchStiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_shearStiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<DistanceConstraint, Real, &DistanceConstraint::m_stiffness>(stiffness);
@@ -234,6 +238,7 @@ void TW_CALL setPoissonRatio(const void* value, void* clientData)
 {
 	poissonRatio = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_poissonRatio>(poissonRatio);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_poissonRatio>(poissonRatio);
 }
 
 void TW_CALL getNormalizeStretch(void* value, void* clientData)

Demos/DistanceFieldDemos/DeformableCollisionDemo.cpp

@@ -88,8 +88,8 @@ int main( int argc, char **argv )
 
 	TwType enumType2 = TwDefineEnum("SimulationMethodType", NULL, 0);
 	TwAddVarCB(MiniGL::getTweakBar(), "SimulationMethod", enumType2, setSimulationMethod, getSimulationMethod, &simulationMethod,
-			" label='Simulation method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics (no inversion handling)}, 4 {Shape matching (no inversion handling)}, 5 {XPBD volume constraints}' group=Simulation");
-	TwAddVarCB(MiniGL::getTweakBar(), "stiffness", TW_TYPE_REAL, setStiffness, getStiffness, model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Solid' ");
+			" label='Simulation method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {FEM based XPBD}, \
+			4 {Strain based dynamics (no inversion handling)}, 5 {Shape matching (no inversion handling)}, 6 {XPBD volume constraints}' group=Simulation");	TwAddVarCB(MiniGL::getTweakBar(), "stiffness", TW_TYPE_REAL, setStiffness, getStiffness, model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Solid' ");
 	TwAddVarCB(MiniGL::getTweakBar(), "poissonRatio", TW_TYPE_REAL, setPoissonRatio, getPoissonRatio, model, " label='Poisson ratio'  min=0.0 step=0.1 precision=4 group='Solid' ");
 	TwAddVarCB(MiniGL::getTweakBar(), "normalizeStretch", TW_TYPE_BOOL32, setNormalizeStretch, getNormalizeStretch, model, " label='Normalize stretch' group='Solid' ");
 	TwAddVarCB(MiniGL::getTweakBar(), "normalizeShear", TW_TYPE_BOOL32, setNormalizeShear, getNormalizeShear, model, " label='Normalize shear' group='Solid' ");
@@ -276,11 +276,13 @@ void createMesh()
 
 	// init constraints
 	stiffness = 1.0;
-	if (simulationMethod == 5)
+	if (simulationMethod == 3)
+		stiffness = 1000000;
+	if (simulationMethod == 6)
 		stiffness = 100000;
 
 	volumeStiffness = 1.0;
-	if (simulationMethod == 5)
+	if (simulationMethod == 6)
 		volumeStiffness = 100000;
 
 	ParticleData& pd = model->getParticles();
@@ -301,6 +303,7 @@ void TW_CALL setStiffness(const void* value, void* clientData)
 {
 	stiffness = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_stiffness>(stiffness);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_stiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_stretchStiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_shearStiffness>(stiffness);
 	((SimulationModel*)clientData)->setConstraintValue<DistanceConstraint, Real, &DistanceConstraint::m_stiffness>(stiffness);
@@ -334,6 +337,7 @@ void TW_CALL setPoissonRatio(const void* value, void* clientData)
 {
 	poissonRatio = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_poissonRatio>(poissonRatio);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_poissonRatio>(poissonRatio);
 }
 
 void TW_CALL getNormalizeStretch(void* value, void* clientData)

Demos/SceneLoaderDemo/SceneLoaderDemo.cpp

@@ -155,7 +155,8 @@ int main( int argc, char **argv )
 	{
 		TwType enumType3 = TwDefineEnum("SolidSimulationMethodType", NULL, 0);
 		TwAddVarCB(MiniGL::getTweakBar(), "SolidSimulationMethod", enumType3, setSolidSimulationMethod, getSolidSimulationMethod, &solidSimulationMethod,
-			" label='Solid sim. method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics (no inversion handling)}, 4 {Shape matching (no inversion handling)}, 5 {XPBD volume constraints}' group=Simulation");
+			" label='Solid sim. method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {FEM based XPBD}, \
+			4 {Strain based dynamics (no inversion handling)}, 5 {Shape matching (no inversion handling)}, 6 {XPBD volume constraints}' group=Simulation");
 		TwAddVarCB(MiniGL::getTweakBar(), "stiffness", TW_TYPE_REAL, setSolidStiffness, getSolidStiffness, model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Solid' ");
 		TwAddVarCB(MiniGL::getTweakBar(), "poissonRatio", TW_TYPE_REAL, setSolidPoissonRatio, getSolidPoissonRatio, model, " label='Poisson ratio'  min=0.0 step=0.1 precision=4 group='Solid' ");
 		TwAddVarCB(MiniGL::getTweakBar(), "normalizeStretch", TW_TYPE_BOOL32, setSolidNormalizeStretch, getSolidNormalizeStretch, model, " label='Normalize stretch' group='Solid' ");
@@ -286,11 +287,13 @@ void initTetModelConstraints()
 	// init constraints
 	SimulationModel *model = Simulation::getCurrent()->getModel();
 	solidStiffness = 1.0;
-	if (solidSimulationMethod == 5)
+	if (solidSimulationMethod == 3)
+		solidStiffness = 1000000;
+	if (solidSimulationMethod == 6)
 		solidStiffness = 100000;
 
 	volumeStiffness = 1.0;
-	if (solidSimulationMethod == 5)
+	if (solidSimulationMethod == 6)
 		volumeStiffness = 100000;
 	for (unsigned int cm = 0; cm < model->getTetModels().size(); cm++)
 	{
@@ -1293,6 +1296,7 @@ void TW_CALL setSolidStiffness(const void* value, void* clientData)
 {
 	solidStiffness = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_stiffness>(solidStiffness);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_stiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_stretchStiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_shearStiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<DistanceConstraint, Real, &DistanceConstraint::m_stiffness>(solidStiffness);
@@ -1326,6 +1330,7 @@ void TW_CALL setSolidPoissonRatio(const void* value, void* clientData)
 {
 	solidPoissonRatio = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_poissonRatio>(solidPoissonRatio);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_poissonRatio>(solidPoissonRatio);
 }
 
 void TW_CALL getSolidNormalizeStretch(void* value, void* clientData)

Demos/StiffRodsDemos/DirectPositionBasedSolverForStiffRodsDemo.cpp

@@ -166,7 +166,8 @@ int main(int argc, char **argv)
 	{
 		TwType enumType3 = TwDefineEnum("SolidSimulationMethodType", NULL, 0);
 		TwAddVarCB(MiniGL::getTweakBar(), "SolidSimulationMethod", enumType3, setSolidSimulationMethod, getSolidSimulationMethod, &solidSimulationMethod,
-			" label='Solid sim. method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {Strain based dynamics (no inversion handling)}, 4 {Shape matching (no inversion handling)}, 5 {XPBD volume constraints}' group=Simulation");
+			" label='Solid sim. method' enum='0 {None}, 1 {Volume constraints}, 2 {FEM based PBD}, 3 {FEM based XPBD}, \
+			4 {Strain based dynamics (no inversion handling)}, 5 {Shape matching (no inversion handling)}, 6 {XPBD volume constraints}' group=Simulation");
 		TwAddVarCB(MiniGL::getTweakBar(), "stiffness", TW_TYPE_REAL, setSolidStiffness, getSolidStiffness, model, " label='Stiffness'  min=0.0 step=0.1 precision=4 group='Solid' ");
 		TwAddVarCB(MiniGL::getTweakBar(), "poissonRatio", TW_TYPE_REAL, setSolidPoissonRatio, getSolidPoissonRatio, model, " label='Poisson ratio'  min=0.0 step=0.1 precision=4 group='Solid' ");
 		TwAddVarCB(MiniGL::getTweakBar(), "normalizeStretch", TW_TYPE_BOOL32, setSolidNormalizeStretch, getSolidNormalizeStretch, model, " label='Normalize stretch' group='Solid' ");
@@ -298,11 +299,13 @@ void initTetModelConstraints()
 	// init constraints
 	SimulationModel* model = Simulation::getCurrent()->getModel();
 	solidStiffness = 1.0;
-	if (solidSimulationMethod == 5)
+	if (solidSimulationMethod == 3)
+		solidStiffness = 1000000;
+	if (solidSimulationMethod == 6)
 		solidStiffness = 100000;
 
 	volumeStiffness = 1.0;
-	if (solidSimulationMethod == 5)
+	if (solidSimulationMethod == 6)
 		volumeStiffness = 100000;
 	for (unsigned int cm = 0; cm < model->getTetModels().size(); cm++)
 	{
@@ -1375,6 +1378,7 @@ void TW_CALL setSolidStiffness(const void* value, void* clientData)
 {
 	solidStiffness = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_stiffness>(solidStiffness);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_stiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_stretchStiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<StrainTetConstraint, Real, &StrainTetConstraint::m_shearStiffness>(solidStiffness);
 	((SimulationModel*)clientData)->setConstraintValue<DistanceConstraint, Real, &DistanceConstraint::m_stiffness>(solidStiffness);
@@ -1408,6 +1412,7 @@ void TW_CALL setSolidPoissonRatio(const void* value, void* clientData)
 {
 	solidPoissonRatio = *(const Real*)(value);
 	((SimulationModel*)clientData)->setConstraintValue<FEMTetConstraint, Real, &FEMTetConstraint::m_poissonRatio>(solidPoissonRatio);
+	((SimulationModel*)clientData)->setConstraintValue<XPBD_FEMTetConstraint, Real, &XPBD_FEMTetConstraint::m_poissonRatio>(solidPoissonRatio);
 }
 
 void TW_CALL getSolidNormalizeStretch(void* value, void* clientData)

PositionBasedDynamics/PositionBasedDynamics.h

@@ -367,7 +367,6 @@ namespace PBD
 
 
 		// -------------- FEM Based PBD  -----------------------------------------------------
-	private:
 		static void computeGradCGreen(
 			Real restVolume, 
 			const Matrix3r &invRestMat, 
@@ -390,7 +389,11 @@ namespace PBD
 
 
 	public:
-		/** Initialize rest configuration infos which are required by the solver step.
+		/** Implementation of the finite element method described in \n\n
+		* Jan Bender, Dan Koschier, Patrick Charrier and Daniel Weber, \n
+		* "Position-Based Simulation of Continuous Materials", \n
+		* Computers & Graphics 44, 2014\n\n
+		* Initialize rest configuration infos which are required by the solver step.
 		* Recomputation is only necessary when rest shape changes.
 		*/
 		static bool init_FEMTriangleConstraint(		
@@ -401,6 +404,12 @@ namespace PBD
 			Matrix2r &invRestMat
 			);
 
+		/** Implementation of the finite element method described in \n\n
+		* Jan Bender, Dan Koschier, Patrick Charrier and Daniel Weber, \n
+		* "Position-Based Simulation of Continuous Materials", \n
+		* Computers & Graphics 44, 2014\n\n
+		* Solve the continuum mechanical constraint defined for a triangle.
+		*/
 		static bool solve_FEMTriangleConstraint(
 			const Vector3r &p0, Real invMass0,
 			const Vector3r &p1, Real invMass1,
@@ -414,7 +423,11 @@ namespace PBD
 			const Real poissonRatioYX,
 			Vector3r &corr0, Vector3r &corr1, Vector3r &corr2);
 
-		/** Initialize rest configuration infos which are required by the solver step.
+		/** Implementation of the finite element method described in \n\n
+		* Jan Bender, Dan Koschier, Patrick Charrier and Daniel Weber, \n
+		* "Position-Based Simulation of Continuous Materials", \n
+		* Computers & Graphics 44, 2014\n\n
+		* Initialize rest configuration infos which are required by the solver step.
 		* Recomputation is only necessary when rest shape changes.
 		*/
 		static bool init_FEMTetraConstraint(			
@@ -426,6 +439,12 @@ namespace PBD
 			Matrix3r &invRestMat
 			);
 
+		/** Implementation of the finite element method described in \n\n
+		* Jan Bender, Dan Koschier, Patrick Charrier and Daniel Weber, \n
+		* "Position-Based Simulation of Continuous Materials", \n
+		* Computers & Graphics 44, 2014\n\n
+		* Solve the continuum mechanical constraint defined for a tetrahedron.
+		*/
 		static bool solve_FEMTetraConstraint(
 			const Vector3r &p0, Real invMass0,
 			const Vector3r &p1, Real invMass1,
@@ -438,6 +457,7 @@ namespace PBD
 			const bool  handleInversion,
 			Vector3r &corr0, Vector3r &corr1, Vector3r &corr2, Vector3r &corr3);
 
+
 		/** Initialize contact between a particle and a tetrahedron and return
 		* info which is required by the solver step.
 		*

PositionBasedDynamics/XPBD.cpp

@@ -1,4 +1,5 @@
 #include "XPBD.h"
+#include "PositionBasedDynamics.h"
 #include "MathFunctions.h"
 #include <cfloat>
 
@@ -211,3 +212,75 @@ bool XPBD::solve_IsometricBendingConstraint(
 	return false;
 }
 
+
+// ----------------------------------------------------------------------------------------------
+bool XPBD::solve_FEMTetraConstraint(
+	const Vector3r& p0, Real invMass0,
+	const Vector3r& p1, Real invMass1,
+	const Vector3r& p2, Real invMass2,
+	const Vector3r& p3, Real invMass3,
+	const Real restVolume,
+	const Matrix3r& invRestMat,
+	const Real youngsModulus,
+	const Real poissonRatio,
+	const bool  handleInversion,
+	const Real dt,
+	Real& multiplier,
+	Vector3r& corr0, Vector3r& corr1, Vector3r& corr2, Vector3r& corr3)
+{
+	corr0.setZero();
+	corr1.setZero();
+	corr2.setZero();
+	corr3.setZero();
+
+	if (youngsModulus <= 0.0)
+		return true;
+
+	if (poissonRatio < 0.0 || poissonRatio > 0.49)
+		return false;
+
+
+	Real C = 0.0;
+	Vector3r gradC[4];
+	Matrix3r epsilon, sigma;
+	Real volume = (p1 - p0).cross(p2 - p0).dot(p3 - p0) / static_cast<Real>(6.0);
+
+	Real mu = 1.0 / static_cast<Real>(2.0) / (static_cast<Real>(1.0) + poissonRatio);
+	Real lambda = 1.0 * poissonRatio / (static_cast<Real>(1.0) + poissonRatio) / (static_cast<Real>(1.0) - static_cast<Real>(2.0) * poissonRatio);
+
+	if (!handleInversion || volume > 0.0)
+	{
+		PositionBasedDynamics::computeGreenStrainAndPiolaStress(p0, p1, p2, p3, invRestMat, restVolume, mu, lambda, epsilon, sigma, C);
+		PositionBasedDynamics::computeGradCGreen(restVolume, invRestMat, sigma, gradC);
+	}
+	else
+	{
+		PositionBasedDynamics::computeGreenStrainAndPiolaStressInversion(p0, p1, p2, p3, invRestMat, restVolume, mu, lambda, epsilon, sigma, C);
+		PositionBasedDynamics::computeGradCGreen(restVolume, invRestMat, sigma, gradC);
+	}
+
+	C = sqrt(2.0 * C);
+
+	Real sum_normGradC =
+		invMass0 * gradC[0].squaredNorm() +
+		invMass1 * gradC[1].squaredNorm() +
+		invMass2 * gradC[2].squaredNorm() +
+		invMass3 * gradC[3].squaredNorm();
+
+	Real alpha = static_cast<Real>(1.0) / (youngsModulus * dt * dt);
+	sum_normGradC += C * C * alpha;
+
+	if (sum_normGradC < eps)
+		return false;
+
+	// compute scaling factor
+	//const Real s = C / sum_normGradC;
+	const Real s = (C * C + C * alpha * multiplier) / sum_normGradC;
+
+	corr0 = -s * invMass0 * gradC[0];
+	corr1 = -s * invMass1 * gradC[1];
+	corr2 = -s * invMass2 * gradC[2];
+	corr3 = -s * invMass3 * gradC[3];
+
+	return true;
+}

PositionBasedDynamics/XPBD.h

@@ -126,6 +126,27 @@ namespace PBD
 			const Real dt,
 			Real& lambda,
 			Vector3r& corr0, Vector3r& corr1, Vector3r& corr2, Vector3r& corr3);
+
+
+		/** XPBD implementation of the finite element method described in \n\n
+		* Jan Bender, Dan Koschier, Patrick Charrier and Daniel Weber, \n
+		* "Position-Based Simulation of Continuous Materials", \n
+		* Computers & Graphics 44, 2014\n\n
+		* Solve the continuum mechanical constraint defined for a tetrahedron.
+		*/
+		static bool solve_FEMTetraConstraint(
+			const Vector3r& p0, Real invMass0,
+			const Vector3r& p1, Real invMass1,
+			const Vector3r& p2, Real invMass2,
+			const Vector3r& p3, Real invMass3,
+			const Real restVolume,
+			const Matrix3r& invRestMat,
+			const Real youngsModulus,
+			const Real poissonRatio,
+			const bool  handleInversion,
+			const Real dt,
+			Real& lambda,
+			Vector3r& corr0, Vector3r& corr1, Vector3r& corr2, Vector3r& corr3);
 	};
 }