- updated to Eigen 3.4.0

- updated to PBD 2.0.1

Author: janbender

CMake/Common.cmake

@@ -9,20 +9,20 @@ set(CMAKE_RELWITHDEBINFO_POSTFIX "_rd")
 set(CMAKE_MINSIZEREL_POSTFIX "_ms")
 
 include(CMakeDependentOption)
+include(avx)
 
 OPTION(USE_DOUBLE_PRECISION "Use double precision"	OFF)
 if (USE_DOUBLE_PRECISION)
 	add_definitions( -DUSE_DOUBLE)	
-	message("AVX is only supported for single precision.")
-endif (USE_DOUBLE_PRECISION)
-
-cmake_dependent_option(USE_AVX "Use AVX" ON "NOT USE_DOUBLE_PRECISION" OFF)
-if (USE_AVX)
-	message(STATUS "If your CPU does not support AVX, turn off the USE_AVX flag.")
-	add_definitions(-DUSE_AVX)
+else()
+	set_avx_flags()
+	if (FOUND_AVX2)
+		message(STATUS "Using AVX2")
+		add_definitions(-DUSE_AVX)
+	endif()
 endif()
 
-cmake_dependent_option(USE_PERFORMANCE_OPTIMIZATION "Optimize performance (higher memory consumption)" ON "USE_AVX" OFF)
+cmake_dependent_option(USE_PERFORMANCE_OPTIMIZATION "Optimize performance (higher memory consumption)" ON "FOUND_AVX2" OFF)
 if (USE_PERFORMANCE_OPTIMIZATION)
 	add_definitions( -DUSE_PERFORMANCE_OPTIMIZATION)	
 endif (USE_PERFORMANCE_OPTIMIZATION)
@@ -70,8 +70,8 @@ if (MSVC)
     set(CMAKE_CXX_FLAGS_RELEASE "/MD /MP /Ox /Ob2 /Oi /Ot /fp:fast /D NDEBUG") 
 	set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "/MD /Zi /MP /Ox /Ob2 /Oi /Ot /fp:fast /D NDEBUG") 
 	set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} /MP") 
-	if (USE_AVX)
-		add_compile_options("/arch:AVX")
+	if (FOUND_AVX2)
+		add_compile_options(${AVX_FLAGS})
 	endif()
 endif (MSVC)
 
@@ -86,9 +86,8 @@ if (UNIX OR MINGW)
 		set(CMAKE_CXX_FLAGS_RELWITHDEBINFO "-O3 -DNDEBUG -march=native")
 	endif ()
 	set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG}")
-	if (USE_AVX)
-		add_compile_options("-mavx")
-		add_compile_options("-mfma")
+	if (FOUND_AVX2)
+		add_compile_options(${AVX_FLAGS})
 	endif()
 endif (UNIX OR MINGW)
 if(MINGW)

CMake/SetUpExternalProjects.cmake

@@ -47,12 +47,12 @@ ExternalProject_Add(
 	Ext_PBD
 	PREFIX "${CMAKE_BINARY_DIR}/extern/PositionBasedDynamics"
 	GIT_REPOSITORY https://github.com/InteractiveComputerGraphics/PositionBasedDynamics.git
-	GIT_TAG "67cea4478d58de55f64f78fcfea16629bfb79152"
+	GIT_TAG "544cc0cba4db8cabe2217cd48b87fb095848ea92"
 	INSTALL_DIR ${ExternalInstallDir}/PositionBasedDynamics
 	DEPENDS Ext_GenericParameters Ext_Discregrid
 	CMAKE_ARGS -DCMAKE_BUILD_TYPE=${EXT_CMAKE_BUILD_TYPE}
 	-DCMAKE_INSTALL_PREFIX:PATH=${ExternalInstallDir}/PositionBasedDynamics
-	-DPBD_NO_DEMOS:BOOL=1
+	-DPBD_LIBS_ONLY:BOOL=1
 	-DPBD_EXTERNALINSTALLDIR:PATH=${ExternalInstallDir}
 	-DUSE_DOUBLE_PRECISION:BOOL=${USE_DOUBLE_PRECISION}
 	-DGenericParameters_INCLUDE_DIR:PATH=${GENERICPARAMETERS_INCLUDE_DIR}

CMakeLists.txt

@@ -31,9 +31,12 @@ add_subdirectory(extern/partio)
 add_subdirectory(extern/md5)
 add_subdirectory(extern/tinyexpr)
 if (NOT SPH_LIBS_ONLY)
-	add_subdirectory(extern/AntTweakBar)
 	add_subdirectory(extern/glfw)
-	add_subdirectory(extern/imgui)
+	if (USE_IMGUI)
+		add_subdirectory(extern/imgui)
+	else()
+		add_subdirectory(extern/AntTweakBar)
+	endif()
 	if (USE_PYTHON_BINDINGS)
 		add_subdirectory(extern/pybind)
 	endif ()

Changelog.txt

@@ -1,3 +1,7 @@
+2.11.3
+	- updated to Eigen 3.4.0
+	- updated to PBD 2.0.1
+	
 2.11.2
 	- fixed command line parameter evaluation
 	- updated documentation

SPlisHSPlasH/Elasticity/Elasticity_Kugelstadt2021.cpp

@@ -629,8 +629,9 @@ void Elasticity_Kugelstadt2021::initFactorization(std::shared_ptr<Factorization>
 	unsigned int row_index = 0;
 	for (int i = 0; i < (int)numParticles; i++)
 	{
-		int particleIndex = group[i];
-		const unsigned int i0 = m_current_to_initial_index[particleIndex];
+		int particleIndex0 = group[i];
+		const unsigned int i0 = particleIndex0;
+		unsigned int particleIndex = m_initial_to_current_index[i0];
 
 		const double restVolumes_id = m_restVolumes[i];
 
@@ -651,7 +652,7 @@ void Elasticity_Kugelstadt2021::initFactorization(std::shared_ptr<Factorization>
 		for (unsigned int j = 0; j < numNeighbors; j++)
 		{
 			const unsigned int neighborIndex0 = m_initialNeighbors[i0][j];
-			const unsigned int neighborIndex = m_initial_to_current_index[m_initialNeighbors[i0][j]];
+			const unsigned int neighborIndex = m_initial_to_current_index[neighborIndex0];
 			const Eigen::Vector3d xj0d = m_model->getPosition0(neighborIndex0).cast<double>();
 
 			const Eigen::Vector3d correctedKernelGradient = m_L[particleIndex].cast<double>() * sim->gradW((xi0d - xj0d).cast<Real>()).cast<double>();
@@ -680,7 +681,7 @@ void Elasticity_Kugelstadt2021::initFactorization(std::shared_ptr<Factorization>
 			for (unsigned int k = 0; k < numNeighbors; k++)
 			{
 				const unsigned int kIndex0 = m_initialNeighbors[i0][k];
-				const unsigned int kIndex = m_initial_to_current_index[m_initialNeighbors[i0][k]];
+				const unsigned int kIndex = m_initial_to_current_index[kIndex0];
 				const Eigen::Vector3d xk0d = m_model->getPosition0(kIndex0).cast<double>();
 
 				const Eigen::Vector3d correctedKernelGradientd = m_L[particleIndex].cast<double>() * sim->gradW((xi0d - xk0d).cast<Real>()).cast<double>();
@@ -810,6 +811,7 @@ void Elasticity_Kugelstadt2021::performNeighborhoodSearchSort()
 	d.sort_field(&m_current_to_initial_index[0]);
 	d.sort_field(&m_L[0]);
 	d.sort_field(&m_restVolumes[0]);
+	d.sort_field(&m_vDiff[0]);
 
 	for (unsigned int i = 0; i < numPart; i++)
 		m_initial_to_current_index[m_current_to_initial_index[i]] = i;
@@ -845,9 +847,9 @@ void Elasticity_Kugelstadt2021::computeMatrixL()
 			//////////////////////////////////////////////////////////////////////////
 			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 unsigned int neighborIndex = m_initial_to_current_index[neighborIndex0];
 
 				const Vector3r &xj0 = m_model->getPosition0(neighborIndex0);
 				const Vector3r xj_xi_0 = xj0 - xi0;
@@ -1002,7 +1004,7 @@ void Elasticity_Kugelstadt2021::rotationMatricesToAVXQuaternions()
 			{
 				const int count = std::min(numParticles - i * 8, 8);
 
-				// store the deformation gradient of 8 particles in avx vectors
+				// store the rotation matrices of 8 particles in avx vectors
 				int idx[8];
 				for (int j = 0; j < count; j++)
 					idx[j] = m_initial_to_current_index[group[8 * i + j]];
@@ -1011,9 +1013,9 @@ void Elasticity_Kugelstadt2021::rotationMatricesToAVXQuaternions()
 
 				Quaternionr q[8];
 				for (auto j = 0; j < count; j++)
-					q[j] = Quaternionr(m_rotations[idx[j]].transpose());
+					q[j] = Quaternionr(m_rotations[idx[j]]).normalized();
 				for (auto j = count; j < 8; j++)
-					q[j] = Quaternionr();
+					q[j] = Quaternionr(1,0,0,0);
 
 				Quaternion8f& q_avx = quats[i];
 				q_avx.set(q);
@@ -1128,9 +1130,9 @@ void Elasticity_Kugelstadt2021::computeRotations()
 				R3.store(r2.data());
 				for (auto j = 0; j < count; j++)
 				{
-					m_rotations[idx[j]].row(0) = r0[j];
-					m_rotations[idx[j]].row(1) = r1[j];
-					m_rotations[idx[j]].row(2) = r2[j];
+					m_rotations[idx[j]].col(0) = r0[j];
+					m_rotations[idx[j]].col(1) = r1[j];
+					m_rotations[idx[j]].col(2) = r2[j];
 				}
 			}
 		}
@@ -1588,9 +1590,9 @@ void Elasticity_Kugelstadt2021::computeRHS(VectorXr & rhs)
 				//const Real trace = RTF(0, 0) + RTF(1, 1) + RTF(2, 2) - 3.0;
 
 				// short form of: trace(R^T F - I) 
-				const Real trace =	m_rotations[i].col(0).dot(m_F[i].col(0)) + 
-									m_rotations[i].col(1).dot(m_F[i].col(1)) + 
-									m_rotations[i].col(2).dot(m_F[i].col(2)) - static_cast<Real>(3.0);
+				const Real trace =	m_rotations[i].row(0).dot(m_F[i].col(0)) + 
+									m_rotations[i].row(1).dot(m_F[i].col(1)) + 
+									m_rotations[i].row(2).dot(m_F[i].col(2)) - static_cast<Real>(3.0);
 
 				m_stress[i] = m_lambda * trace;
 			}
@@ -2215,7 +2217,7 @@ void Elasticity_Kugelstadt2021::precomputeValues()
 		#pragma omp for schedule(static) 
 		for (int i = 0; i < (int)numParticles; i++)
 		{
-			m_RL[i] = m_rotations[i] * m_L[i];
+			m_RL[i] = m_rotations[i].transpose() * m_L[i];
 		}
 		#pragma omp for schedule(static) 
 		for (int i = 0; i < (int)numParticles; i++)

SPlisHSPlasH/Elasticity/Elasticity_Peer2018.h

@@ -10,10 +10,10 @@ namespace SPH
 {
 	/** \brief This class implements the implicit SPH formulation for 
 	* incompressible linearly elastic solids introduced
-	* by Peer et al. [PGBT17].
+	* by Peer et al. [PGBT18].
 	*
 	* References:
-	* - [PGBT17] Andreas Peer, Christoph Gissler, Stefan Band, and Matthias Teschner. An implicit SPH formulation for incompressible linearly elastic solids. Computer Graphics Forum, 2017. URL: http://dx.doi.org/10.1111/cgf.13317
+	* - [PGBT18] Andreas Peer, Christoph Gissler, Stefan Band, and Matthias Teschner. An implicit SPH formulation for incompressible linearly elastic solids. Computer Graphics Forum, 2018. URL: http://dx.doi.org/10.1111/cgf.13317
 	*/
 	class Elasticity_Peer2018 : public ElasticityBase
 	{

SPlisHSPlasH/FluidModel.cpp

@@ -493,11 +493,13 @@ void FluidModel::setSurfaceTensionMethod(const unsigned int val)
 	m_surfaceTension = nullptr;
 
 	m_surfaceTensionMethod = stm;
-	int method_id = getValue<int>(FluidModel::SURFACE_TENSION_METHOD);
 	for (unsigned int i = 0; i < stMethods.size(); i++)
 	{
-		if (stMethods[i].m_id == method_id)
+		if (stMethods[i].m_id == m_surfaceTensionMethod)
+		{
 			m_surfaceTension = static_cast<SurfaceTensionBase*>(stMethods[i].m_creator(this));
+			break;
+		}
 	}
 
 	if (m_surfaceTension != nullptr)
@@ -585,10 +587,9 @@ void FluidModel::setVorticityMethod(const unsigned int val)
 
 	m_vorticityMethod = vm;
 
-	int method_id = getValue<int>(FluidModel::VORTICITY_METHOD);
 	for (unsigned int i = 0; i < vorticityMethods.size(); i++)
 	{
-		if (vorticityMethods[i].m_id == method_id)
+		if (vorticityMethods[i].m_id == m_vorticityMethod)
 			m_vorticity = static_cast<VorticityBase*>(vorticityMethods[i].m_creator(this));
 	}
 
@@ -630,10 +631,9 @@ void FluidModel::setDragMethod(const unsigned int val)
 
 	m_dragMethod = dm;
 
-	int method_id = getValue<int>(FluidModel::DRAG_METHOD);
 	for (unsigned int i = 0; i < dragMethods.size(); i++)
 	{
-		if (dragMethods[i].m_id == method_id)
+		if (dragMethods[i].m_id == m_dragMethod)
 			m_drag = static_cast<DragBase*>(dragMethods[i].m_creator(this));
 	}
 
@@ -674,11 +674,9 @@ void FluidModel::setElasticityMethod(const unsigned int val)
 	m_elasticity = nullptr;
 
 	m_elasticityMethod = em;
-
-	int method_id = getValue<int>(FluidModel::ELASTICITY_METHOD);
 	for (unsigned int i = 0; i < elasticityMethods.size(); i++)
 	{
-		if (elasticityMethods[i].m_id == method_id)
+		if (elasticityMethods[i].m_id == m_elasticityMethod)
 			m_elasticity = static_cast<ElasticityBase*>(elasticityMethods[i].m_creator(this));
 	}
 

SPlisHSPlasH/RigidBodyObject.h

@@ -19,7 +19,7 @@ namespace SPH
 
 		virtual bool isDynamic() const = 0;
 		bool isAnimated() const { return m_isAnimated; }
-		virtual void setIsAnimated(const bool b) { m_isAnimated = true; }
+		virtual void setIsAnimated(const bool b) { m_isAnimated = b; }
 
 		virtual Real const getMass() const = 0;
 		virtual Vector3r const& getPosition() const = 0;

SPlisHSPlasH/TimeStep.cpp

@@ -525,7 +525,6 @@ void TimeStep::precomputeValues()
 		{
 			for (unsigned int pid = 0; pid < nFluids; pid++)
 			{
-				FluidModel* fm_neighbor = sim->getFluidModelFromPointSet(pid);
 				const unsigned int maxN = sim->numberOfNeighbors(fluidModelIndex, pid, i);
 
 				// same phase
@@ -551,7 +550,7 @@ void TimeStep::precomputeValues()
 			{
 				const Vector3r& xi = model->getPosition(i);
 				const Vector3f8 xi_avx(xi);
-				unsigned int base = precomputed_indices[i];
+				const unsigned int base = precomputed_indices[i];
 				unsigned int idx = 0;
 				forall_fluid_neighbors_avx(
 					const Scalarf8 Vj_avx = convert_zero(fm_neighbor->getVolume(0), count);

SPlisHSPlasH/Utilities/AVX_math.h

@@ -957,6 +957,23 @@ class Quaternion8f
 		}
 	}
 
+	inline void store(Quaternionr *qf) const
+	{
+		float x[8], y[8], z[8], w[8];
+		q[0].store(x);
+		q[1].store(y);
+		q[2].store(z);
+		q[3].store(w);
+
+		for (int i = 0; i < 8; i++)
+		{
+			qf[i].x() = x[i];
+			qf[i].y() = y[i];
+			qf[i].z() = z[i];
+			qf[i].w() = w[i];
+		}
+	}
+
 	inline void set(const Quaternionr *qf)
 	{
 		float x[8], y[8], z[8], w[8];