Merge branch 'InteractiveComputerGraphics:master' into fix-macos

Author: ruberith

CMake/Common.cmake

@@ -10,6 +10,7 @@ set(CMAKE_MINSIZEREL_POSTFIX "_ms")
 
 include(CMakeDependentOption)
 
+
 OPTION(USE_DOUBLE_PRECISION "Use double precision"	OFF)
 if (USE_DOUBLE_PRECISION)
 	add_definitions( -DUSE_DOUBLE)	
@@ -18,11 +19,28 @@ 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 or Accelerate, turn off the USE_AVX flag.")
-	add_definitions(-DUSE_AVX)
+	set(FOUND_ACCELERATE OFF)
+	if (APPLE)
+		find_library(ACCELERATE Accelerate)
+	endif()
+	if (ACCELERATE)
+		set(FOUND_ACCELERATE ON)
+		message(STATUS "Using Accelerate")
+		link_libraries(${ACCELERATE})
+	else()
+		include(avx)
+		set_avx_flags()
+		if (FOUND_AVX2)
+			message(STATUS "Using AVX2")
+		endif()
+	endif()
+	if (FOUND_ACCELERATE OR FOUND_AVX2)
+		set(FOUND_SIMD ON)
+		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_SIMD" OFF)
 if (USE_PERFORMANCE_OPTIMIZATION)
 	add_definitions( -DUSE_PERFORMANCE_OPTIMIZATION)	
 endif (USE_PERFORMANCE_OPTIMIZATION)
@@ -74,8 +92,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)
 
@@ -93,13 +111,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)
-		if (APPLE)
-			link_libraries("-framework Accelerate")
-		else()
-			add_compile_options("-mavx")
-			add_compile_options("-mfma")
-		endif()
+	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/ruberith/PositionBasedDynamics.git
-	GIT_TAG "66081a7fe5b62a7e44553d49a7c3f12b38275544"
+	GIT_TAG "3bdc9079b8354a8ee2c25c3d842f4e8c96077eb3"
 	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}

CMake/avx.cmake

@@ -0,0 +1,71 @@
+macro(set_avx_flags)
+	set(AVX_FLAGS)
+
+	include(CheckCXXSourceRuns)
+	set(CMAKE_REQUIRED_FLAGS)
+
+	# AVX
+	if(MSVC AND NOT MSVC_VERSION LESS 1600)
+		set(CMAKE_REQUIRED_FLAGS "/arch:AVX")
+	elseif(UNIX OR MINGW)
+		set(CMAKE_REQUIRED_FLAGS "-mavx")
+	endif()
+
+	check_cxx_source_runs("
+		#include <immintrin.h>
+		int main()
+		{
+			float v[8] = { 1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f };
+			float r[8];
+			__m256 first = _mm256_loadu_ps(v);
+			__m256 second = _mm256_setr_ps(0.5f, 0.5f, 0.5f, 0.5f, 0.5f, 0.5f, 0.5f, 0.5f);
+			__m256 result = _mm256_add_ps(first, second);
+			_mm256_storeu_ps(r, result);
+
+			for( int i = 0; i < 8; i++)
+				if ((i+1)*1.0f+0.5f != r[i])
+					return -1;
+			return 0;
+		}"
+			FOUND_AVX)
+
+	# AVX2
+	if(MSVC AND NOT MSVC_VERSION LESS 1800)
+		set(CMAKE_REQUIRED_FLAGS "/arch:AVX2")
+	elseif(UNIX OR MINGW)
+		set(CMAKE_REQUIRED_FLAGS "-mavx2")
+	endif()
+
+	check_cxx_source_runs("
+		#include <immintrin.h>
+		int main()
+		{
+			int v[8] = {10, 20, 30, 40, 50, 60, 70, 80};
+			int r[8];
+			__m256i first = _mm256_loadu_si256((__m256i*)v);
+			__m256i second = _mm256_set_epi32(5, 5, 5, 5, 5, 5, 5, 5);
+			__m256i result = _mm256_add_epi32(first, second);
+			_mm256_storeu_si256((__m256i*)r, result);
+
+			for( int i = 0; i < 8; i++)
+				if ((i+1)*10+5 != r[i])
+					return -1;
+			return 0;
+		}"
+			FOUND_AVX2)
+
+	# set compiler flags
+	if (FOUND_AVX2)
+		if(MSVC)
+			set(AVX_FLAGS "/arch:AVX2")
+		elseif(UNIX OR MINGW)
+			set(AVX_FLAGS "-mavx2;-mfma")
+		endif()
+	elseif (FOUND_AVX)
+		if(MSVC)
+			set(AVX_FLAGS "/arch:AVX")
+		elseif(UNIX OR MINGW)
+			set(AVX_FLAGS "-mavx")
+		endif()
+	endif()
+endmacro()

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