Official implementation of the TOG paper "Implicit Surface Tension fo…

…r SPH Fluid Simulation" (#291)

* start porting surface tension

add scenes

* add python bindings

* bump version and fix other surface tension scenes

* fix non-avx surface tension model

* bump to next major version

SPlisHSPlasH/CMakeLists.txt

@@ -75,13 +75,15 @@ set(SURFACETENSION_HEADER_FILES
 	SurfaceTension/SurfaceTension_Becker2007.h
 	SurfaceTension/SurfaceTension_Akinci2013.h
 	SurfaceTension/SurfaceTension_He2014.h
+	SurfaceTension/SurfaceTension_Jeske2023.h
 	)
 
 set(SURFACETENSION_SOURCE_FILES
 	SurfaceTension/SurfaceTensionBase.cpp
 	SurfaceTension/SurfaceTension_Becker2007.cpp
 	SurfaceTension/SurfaceTension_Akinci2013.cpp
 	SurfaceTension/SurfaceTension_He2014.cpp
+	SurfaceTension/SurfaceTension_Jeske2023.cpp
 	)	
 
 set(VISCOSITY_HEADER_FILES

SPlisHSPlasH/NonPressureForceRegistration.cpp

@@ -20,6 +20,7 @@
 #include "SurfaceTension/SurfaceTension_Becker2007.h"
 #include "SurfaceTension/SurfaceTension_Akinci2013.h"
 #include "SurfaceTension/SurfaceTension_He2014.h"
+#include "SurfaceTension/SurfaceTension_Jeske2023.h"
 #ifdef USE_THIRD_PARTY_METHODS
 #include "SurfaceTension/SurfaceTension_ZorillaRitter2020.h"
 #endif
@@ -42,6 +43,7 @@ void Simulation::registerNonpressureForces()
 	addSurfaceTensionMethod("Becker & Teschner 2007", SurfaceTension_Becker2007::creator);
 	addSurfaceTensionMethod("Akinci et al. 2013", SurfaceTension_Akinci2013::creator);
 	addSurfaceTensionMethod("He et al. 2014", SurfaceTension_He2014::creator);
+	addSurfaceTensionMethod("Jeske et al. 2023", SurfaceTension_Jeske2023::creator);
 #ifdef USE_THIRD_PARTY_METHODS
 	addSurfaceTensionMethod("Zorilla, Ritter, et al. 2020", SurfaceTension_ZorillaRitter2020::creator);
 #endif

SPlisHSPlasH/SurfaceTension/SurfaceTension_Jeske2023.h

@@ -0,0 +1,115 @@
+#ifndef __SurfaceTension_Jeske2023_Surface_Tension_h__
+#define __SurfaceTension_Jeske2023_Surface_Tension_h__
+
+#include "SPlisHSPlasH/Common.h"
+#include "SPlisHSPlasH/FluidModel.h"
+#include "SurfaceTensionBase.h"
+#include "SPlisHSPlasH/Utilities/MatrixFreeSolver.h"
+#include "Utilities/Logger.h"
+
+namespace SPH
+{
+	/** \brief This class implements the implicit surface tension method 
+	* by Jeske et al. 2023 [JWL+23].
+	*
+	* References:
+	* - [JWL+23] Jeske, Stefan Rhys, Lukas Westhofen, Fabian Löschner, José Antonio Fernández-Fernández, and Jan Bender. “Implicit Surface Tension for SPH Fluid Simulation.” ACM Transactions on Graphics, November 7, 2023. URL: https://doi.org/10.1145/3631936.
+
+	*/
+	class SurfaceTension_Jeske2023 : public SurfaceTensionBase
+	{
+	protected:
+        Real m_viscosity;
+        Real m_boundaryViscosity;
+
+		unsigned int m_maxIter;
+		Real m_maxError;
+		unsigned int m_iterations;
+		std::vector<Vector3r> m_vDiff;
+		std::vector<Real> m_gradRho;
+        std::vector<Real> m_surfaceEnergy;
+        std::vector<Real> m_color;
+        std::vector<Vector3r> m_colorGrad;
+        std::vector<Vector3r> m_nonlinearAcc;
+        std::vector<Vector3r> m_nonlinearRes;
+        std::vector<Vector3r> m_nonlinearGrad;
+
+		Real m_tangentialDistanceFactor;
+		bool m_weakPhaseCoupling;
+        Real m_xsph;
+
+		typedef Eigen::ConjugateGradient<MatrixReplacement, Eigen::Lower | Eigen::Upper, Eigen::IdentityPreconditioner> Solver;
+
+		Solver m_solver;
+
+
+		virtual void initParameters();
+
+	public:
+		static int ITERATIONS;
+		static int MAX_ITERATIONS;
+		static int MAX_ERROR;
+        static int VISCOSITY_COEFFICIENT;
+		static int VISCOSITY_COEFFICIENT_BOUNDARY;
+        static int XSPH;
+
+		SurfaceTension_Jeske2023(FluidModel *model);
+		virtual ~SurfaceTension_Jeske2023(void);
+
+		static NonPressureForceBase* creator(FluidModel* model) { return new SurfaceTension_Jeske2023(model); }
+
+		virtual void step();
+		virtual void reset();
+
+		virtual void performNeighborhoodSearchSort();
+
+		static void matrixVecProd(const Real* vec, Real *result, void *userData);
+
+		FORCE_INLINE const Vector3r& getVDiff(const unsigned int i) const
+		{
+			return m_vDiff[i];
+		}
+
+		FORCE_INLINE Vector3r& getVDiff(const unsigned int i)
+		{
+			return m_vDiff[i];
+		}
+
+		FORCE_INLINE void setVDiff(const unsigned int i, const Vector3r& val)
+		{
+			m_vDiff[i] = val;
+		}
+
+        FORCE_INLINE const Real& getDensityGrad(const unsigned int i) const
+        {
+            return m_gradRho[i];
+        }
+
+        FORCE_INLINE Real& getDensityGrad(const unsigned int i)
+        {
+            return m_gradRho[i];
+        }
+
+        FORCE_INLINE void setDensityGrad(const unsigned int i, const Real& val)
+        {
+            m_gradRho[i] = val;
+        }
+
+        void computeRHS(VectorXr &b, VectorXr &g);
+
+        void applyForces(const VectorXr &x);
+
+        Real getMaxSolverError(){ return m_maxError; }
+        void setMaxSolverError(Real error){ m_maxError = error; }
+
+        bool getWeakCoupling(){ return m_weakPhaseCoupling; }
+        void setWeakCoupling(bool val){ m_weakPhaseCoupling = val; }
+
+        bool getViscosity(){ return m_viscosity; }
+        void setViscosity(Real val){ m_viscosity = val; }
+
+		void computeDensityGradient();
+    };
+}
+
+#endif

SPlisHSPlasH/Utilities/AVX_math.h

@@ -153,6 +153,10 @@ static inline Scalarf8 operator >= (Scalarf8 const & a, Scalarf8 const & b) {
 	return _mm256_cmp_ps(b.v, a.v, 2);
 }
 
+static inline Scalarf8 min(Scalarf8 const & a, Scalarf8 const & b) {
+	return _mm256_min_ps(a.v, b.v);
+}
+
 static inline Scalarf8 max(Scalarf8 const & a, Scalarf8 const & b) {
 	return _mm256_max_ps(a.v, b.v);
 }

data/Scenes/SurfaceTension_BoxHole_ZR2020.json

@@ -39,7 +39,7 @@
 		{
 			"id": "Fluid", 
 			"surfaceTension": 0.15,
-			"surfaceTensionMethod": 4,
+			"surfaceTensionMethod": 5,
 			"viscosity": 0.01,
 			"viscosityMethod": 1,
 			"surfTZRversion": 1,

data/Scenes/SurfaceTension_BreakDamZR2020.json

@@ -24,7 +24,7 @@
 		{
 			"id": "Fluid",
 			"surfaceTension": 0.1,
-			"surfaceTensionMethod": 4,
+			"surfaceTensionMethod": 5,
 			"viscosity": 0.01,
 			"viscosityMethod": 1,
 			"surfTZRversion": 0,