Oif corrections (#3385)

This ports back fixes from the object in fluid development branch. In particular, the bending force between two triangeles becomes torque-free with this.
The test in this PR fails without the correctiosn in oif_local_force.hpp

maintainer/configs/maxset.hpp

@@ -75,7 +75,5 @@
 
 #define VIRTUAL_SITES_RELATIVE
 #define VIRTUAL_SITES_INERTIALESS_TRACERS
-#define OIF_GLOBAL_FORCES
-#define OIF_LOCAL_FORCES
 
 #define ADDITIONAL_CHECKS

samples/object_in_fluid/motivation.py

@@ -22,7 +22,7 @@
 import espressomd
 
 required_features = ["LB_BOUNDARIES", "EXTERNAL_FORCES", "SOFT_SPHERE",
-                     "OIF_GLOBAL_FORCES", "OIF_LOCAL_FORCES", "MASS"]
+                     "MASS"]
 espressomd.assert_features(required_features)
 
 from espressomd import lbboundaries

src/config/features.def

@@ -91,8 +91,6 @@ OVERLAPPED
 THOLE                           requires ELECTROSTATICS
 
 /* Fluid-Structure Interactions (object in fluid) */
-OIF_LOCAL_FORCES
-OIF_GLOBAL_FORCES
 
 /* Immersed-Boundary Bayreuth version */
 SCAFACOS_DIPOLES                requires SCAFACOS

src/core/forces.cpp

@@ -130,7 +130,6 @@ void force_calc(CellStructure &cell_structure) {
 
   Constraints::constraints.add_forces(particles, sim_time);
 
-#ifdef OIF_GLOBAL_FORCES
   if (max_oif_objects) {
     double area_volume[2]; // There are two global quantities that need to be
     // evaluated: object's surface and object's volume. One
@@ -144,7 +143,6 @@ void force_calc(CellStructure &cell_structure) {
       add_oif_global_forces(area_volume, i, particles);
     }
   }
-#endif
 
   // Must be done here. Forces need to be ghost-communicated
   immersed_boundaries.volume_conservation();

src/core/forces_inline.hpp

@@ -489,11 +489,9 @@ inline void add_bonded_force(Particle *const p1) {
             angle_cossquare_force(p1->r.p, p2->r.p, p3->r.p, iaparams);
         bond_broken = false;
         break;
-#ifdef OIF_GLOBAL_FORCES
       case BONDED_IA_OIF_GLOBAL_FORCES:
         bond_broken = false;
         break;
-#endif
       case BONDED_IA_TABULATED_ANGLE:
         std::tie(force1, force2, force3) =
             tab_angle_force(p1->r.p, p2->r.p, p3->r.p, iaparams);
@@ -516,14 +514,12 @@ inline void add_bonded_force(Particle *const p1) {
     } // 2 partners (angle bonds...)
     else if (n_partners == 3) {
       switch (type) {
-#ifdef OIF_LOCAL_FORCES
       case BONDED_IA_OIF_LOCAL_FORCES:
         // in OIF nomenclature, particles p2 and p3 are common to both triangles
         std::tie(force1, force2, force3, force4) =
             calc_oif_local(*p1, *p2, *p3, *p4, iaparams);
         bond_broken = false;
         break;
-#endif
       // IMMERSED_BOUNDARY
       case BONDED_IA_IBM_TRIBEND:
         std::tie(force1, force2, force3, force4) =
@@ -605,9 +601,7 @@ inline void add_bonded_force(Particle *const p1) {
       case BONDED_IA_TABULATED_DIHEDRAL:
       case BONDED_IA_DIHEDRAL:
       case BONDED_IA_IBM_TRIBEND:
-#ifdef OIF_LOCAL_FORCES
       case BONDED_IA_OIF_LOCAL_FORCES:
-#endif
         p1->f.f += force1;
         p2->f.f += force2;
         p3->f.f += force3;

src/core/global.cpp

@@ -159,10 +159,8 @@ const std::unordered_map<int, Datafield> fields{
      {langevin.gamma_rotation.data(), Datafield::Type::DOUBLE, 3,
       "langevin.gamma_rotation"}}, /* 55 from thermostat.cpp */
 #endif
-#ifdef OIF_GLOBAL_FORCES
     {FIELD_MAX_OIF_OBJECTS,
      {&max_oif_objects, Datafield::Type::INT, 1, "max_oif_objects"}},
-#endif
     {FIELD_THERMALIZEDBONDS,
      {&n_thermalized_bonds, Datafield::Type::INT, 1,
       "n_thermalized_bonds"}}, /* 56 from thermalized_bond.cpp */

src/core/object-in-fluid/oif_local_forces.hpp

@@ -33,7 +33,7 @@
 #include <utils/Vector.hpp>
 #include <utils/math/triangle_functions.hpp>
 
-/** Set parameters for OIF local forces */
+// set parameters for local forces
 int oif_local_forces_set_params(int bond_type, double r0, double ks,
                                 double kslin, double phi0, double kb,
                                 double A01, double A02, double kal,
@@ -134,20 +134,34 @@ calc_oif_local(Particle const &p2, Particle const &p1, Particle const &p3,
      force for triangle p2,p3,p4 p1 += forceT1; p2 -= 0.5*forceT1+0.5*forceT2;
      p3 -= 0.5*forceT1+0.5*forceT2; p4 += forceT2; */
   if (iaparams.p.oif_local_forces.kb > TINY_OIF_ELASTICITY_COEFFICIENT) {
-    auto const n1 = Utils::get_n_triangle(fp2, fp1, fp3).normalize();
-    auto const n2 = Utils::get_n_triangle(fp2, fp3, fp4).normalize();
-
     auto const phi = Utils::angle_btw_triangles(fp1, fp2, fp3, fp4);
     auto const aa = (phi - iaparams.p.oif_local_forces
                                .phi0); // no renormalization by phi0, to be
                                        // consistent with Krueger and Fedosov
     auto const fac = iaparams.p.oif_local_forces.kb * aa;
-    auto const f = 0.5 * fac * n1 + 0.5 * fac * n2;
 
-    force1 += fac * n1;
-    force2 -= f;
-    force3 -= f;
-    force4 += fac * n2;
+    auto const Nc = Utils::get_n_triangle(
+        fp1, fp2,
+        fp3); // returns (fp2 - fp1)x(fp3 - fp1), thus Nc = (A - C)x(B - C)
+    auto const Nd = Utils::get_n_triangle(
+        fp4, fp3,
+        fp2); // returns (fp3 - fp4)x(fp2 - fp4), thus Nd = (B - D)x(A - D)
+
+    auto const BminA = fp3 - fp2;
+
+    auto const factorFaNc =
+        (fp2 - fp3) * (fp1 - fp3) / BminA.norm() / Nc.norm2();
+    auto const factorFaNd =
+        (fp2 - fp3) * (fp4 - fp3) / BminA.norm() / Nd.norm2();
+    auto const factorFbNc =
+        (fp2 - fp3) * (fp2 - fp1) / BminA.norm() / Nc.norm2();
+    auto const factorFbNd =
+        (fp2 - fp3) * (fp2 - fp4) / BminA.norm() / Nd.norm2();
+
+    force1 -= fac * BminA.norm() / Nc.norm2() * Nc;      // Fc
+    force2 += fac * (factorFaNc * Nc + factorFaNd * Nd); // Fa
+    force3 += fac * (factorFbNc * Nc + factorFbNd * Nd); // Fb
+    force4 -= fac * BminA.norm() / Nd.norm2() * Nd;      // Fc
   }
 
   /* local area

src/python/espressomd/system.pyx

@@ -58,15 +58,12 @@ IF COLLISION_DETECTION == 1:
 
 
 setable_properties = ["box_l", "min_global_cut", "periodicity", "time",
-                      "time_step", "timings", "force_cap"]
+                      "time_step", "timings", "force_cap", "max_oif_objects"]
 
 if VIRTUAL_SITES:
     setable_properties.append("_active_virtual_sites_handle")
 
 
-if OIF_GLOBAL_FORCES:
-    setable_properties.append("max_oif_objects")
-
 cdef bool _system_created = False
 
 cdef class System:
@@ -302,19 +299,18 @@ cdef class System:
             def __get__(self):
                 return self._active_virtual_sites_handle.implementation
 
-    IF OIF_GLOBAL_FORCES:
-        property max_oif_objects:
-            """Maximum number of objects as per the object_in_fluid method.
+    property max_oif_objects:
+        """Maximum number of objects as per the object_in_fluid method.
 
-            """
+        """
 
-            def __get__(self):
-                return max_oif_objects
+        def __get__(self):
+            return max_oif_objects
 
-            def __set__(self, v):
-                global max_oif_objects
-                max_oif_objects = v
-                mpi_bcast_parameter(FIELD_MAX_OIF_OBJECTS)
+        def __set__(self, v):
+            global max_oif_objects
+            max_oif_objects = v
+            mpi_bcast_parameter(FIELD_MAX_OIF_OBJECTS)
 
     def change_volume_and_rescale_particles(self, d_new, dir="xyz"):
         """Change box size and rescale particle coordinates.

src/python/object_in_fluid/oif_utils.py

@@ -112,22 +112,12 @@ def angle_btw_triangles(P1, P2, P3, P4):
     """
     n1 = get_triangle_normal(P2, P1, P3)
     n2 = get_triangle_normal(P2, P3, P4)
-    tmp11 = np.dot(n1, n2)
-    tmp11 = tmp11 * abs(tmp11)
-
-    tmp22 = np.dot(n1, n1)
-    tmp33 = np.dot(n2, n2)
-    tmp11 /= (tmp22 * tmp33)
-
-    if tmp11 > 0:
-        tmp11 = np.sqrt(tmp11)
-    else:
-        tmp11 = - np.sqrt(- tmp11)
+    tmp11 = np.dot(n1, n2) / (np.linalg.norm(n1) * np.linalg.norm(n2))
 
     if tmp11 >= 1.0:
-        tmp11 = 0.0
-    elif tmp11 <= -1.:
-        tmp11 = np.pi
+        tmp11 = 1.0
+    elif tmp11 <= -1.0:
+        tmp11 = -1.0
 
     phi = np.pi - math.acos(tmp11)
 

src/utils/include/utils/math/triangle_functions.hpp

@@ -21,6 +21,8 @@
 
 #include "utils/Vector.hpp"
 
+#include <boost/algorithm/clamp.hpp>
+
 #include <cmath>
 
 namespace Utils {
@@ -48,54 +50,46 @@ inline double area_triangle(const Vector3d &P1, const Vector3d &P2,
   return 0.5 * get_n_triangle(P1, P2, P3).norm();
 }
 
-/** This function returns the angle between the triangle p1,p2,p3 and p2,p3,p4.
- *  Be careful, the angle depends on the orientation of the triangles! You need
- *  to be sure that the orientation (direction of normal vector) of p1p2p3
- *  is given by the cross product p2p1 x p2p3. The orientation of p2p3p4 must
- *  be given by p2p3 x p2p4.
- *
- *  Example: p1 = (0,0,1), p2 = (0,0,0), p3=(1,0,0), p4=(0,1,0). The
- *  orientation of p1p2p3 should be in the direction (0,1,0) and indeed: p2p1 x
- *  p2p3 = (0,0,1)x(1,0,0) = (0,1,0) This function is called in the beginning
- *  of the simulation when creating bonds depending on the angle between the
- *  triangles, the bending_force. Here, we determine the orientations by
- *  looping over the triangles and checking the correct orientation. So if you
- *  have the access to the order of particles, you are safe to call this
- *  function with exactly this order. Otherwise you need to check the
- *  orientations.
- */
-template <typename T1, typename T2, typename T3, typename T4>
-double angle_btw_triangles(const T1 &P1, const T2 &P2, const T3 &P3,
-                           const T4 &P4) {
-  auto const normal1 = get_n_triangle(P2, P1, P3);
-  auto const normal2 = get_n_triangle(P2, P3, P4);
+/** @brief Returns the angle between two triangles in 3D space
+
 
-  double tmp11;
-  // Now we compute the scalar product of n1 and n2 divided by the norms of n1
-  // and n2
-  // tmp11 = dot(3,normal1,normal2);         // tmp11 = n1.n2
-  tmp11 = normal1 * normal2; // tmp11 = n1.n2
+Returns the angle between two triangles in 3D space given by points P1P2P3 and
+P2P3P4. Note, that the common edge is given as the second and the third
+argument. Here, the angle can have values from 0 to 2 * PI, depending on the
+orientation of the two triangles. So the angle can be convex or concave. For
+each triangle, an inward direction has been defined as the direction of one of
+the two normal vectors. Particularly, for triangle P1P2P3 it is the vector N1 =
+P2P1 x P2P3 and for triangle P2P3P4 it is N2 = P2P3 x P2P4. The method first
+computes the angle between N1 and N2, which gives always value between 0 and PI
+and then it checks whether this value must be corrected to a value between PI
+and 2 * PI.
 
-  tmp11 *= fabs(tmp11);                         // tmp11 = (n1.n2)^2
-  tmp11 /= (normal1.norm2() * normal2.norm2()); // tmp11 = (n1.n2/(|n1||n2|))^2
-  if (tmp11 > 0) {
-    tmp11 = std::sqrt(tmp11);
-  } else {
-    tmp11 = -std::sqrt(-tmp11);
-  }
+As an example, consider 4 points A,B,C,D in space given by coordinates A =
+[1,1,1], B = [2,1,1], C = [1,2,1], D = [1,1,2]. We want to determine the angle
+between triangles ABC and ACD. In case that the orientations of the triangle ABC
+is [0,0,1] and orientation of ACD is [1,0,0], then the resulting angle must be
+PI/2.0. To get correct result, note that the common edge is AC, and one must
+call the method as angle_btw_triangles(B,A,C,D). With this call we have ensured
+that N1 = AB x AC (which coincides with [0,0,1]) and N2 = AC x AD (which
+coincides with [1,0,0]). Alternatively, if the orientations of the two triangles
+were the oppisite, the correct call would be angle_btw_triangles(B,C,A,D) so
+that N1 = CB x CA and N2 = CA x CD.
 
-  if (tmp11 >= 1.) {
-    tmp11 = 0.0;
-  } else if (tmp11 <= -1.) {
-    tmp11 = M_PI;
-  }
-  auto const phi =
-      M_PI - std::acos(tmp11); // The angle between the faces (not considering
-                               // the orientation, always less or equal to Pi)
-                               // is equal to Pi minus angle between the normals
+*/
+inline double angle_btw_triangles(const Vector3d &P1, const Vector3d &P2,
+                                  const Vector3d &P3, const Vector3d &P4) {
+  auto const normal1 = get_n_triangle(P2, P1, P3);
+  auto const normal2 = get_n_triangle(P2, P3, P4);
+  auto const cosine = boost::algorithm::clamp(
+      normal1 * normal2 / std::sqrt(normal1.norm2() * normal2.norm2()), -1.0,
+      1.0);
+  // The angle between the faces (not considering
+  // the orientation, always less or equal to Pi)
+  // is equal to Pi minus angle between the normals
+  auto const phi = M_PI - std::acos(cosine);
 
   // Now we need to determine, if the angle between two triangles is less than
-  // Pi or more than Pi. To do this we check,
+  // Pi or more than Pi. To do this, we check
   // if the point P4 lies in the halfspace given by triangle P1P2P3 and the
   // normal to this triangle. If yes, we have
   // angle less than Pi, if not, we have angle more than Pi.
@@ -104,8 +98,7 @@ double angle_btw_triangles(const T1 &P1, const T2 &P2, const T3 &P3,
   // Point P1 lies in the plane, therefore d = -(n_x*P1_x + n_y*P1_y + n_z*P1_z)
   // Point P4 lies in the halfspace given by normal iff n_x*P4_x + n_y*P4_y +
   // n_z*P4_z + d >= 0
-  tmp11 = -(normal1[0] * P1[0] + normal1[1] * P1[1] + normal1[2] * P1[2]);
-  if (normal1[0] * P4[0] + normal1[1] * P4[1] + normal1[2] * P4[2] + tmp11 < 0)
+  if (normal1 * P4 - normal1 * P1 < 0)
     return 2 * M_PI - phi;
 
   return phi;

src/utils/tests/triangle_functions_test.cpp

@@ -66,3 +66,25 @@ BOOST_AUTO_TEST_CASE(normal_) {
     BOOST_CHECK_SMALL(normal * P1P3, epsilon);
   }
 }
+
+BOOST_AUTO_TEST_CASE(angle_triangles) {
+  /* Notes from @icimrak from Github #3385
+  As an example, consider 4 points A,B,C,D in space given by coordinates A =
+  [1,1,1], B = [2,1,1], C = [1,2,1], D = [1,1,2]. We want to determine the angle
+  between triangles ABC and ACD. In case that the orientations of the triangle
+  ABC is [0,0,1] and orientation of ACD is [1,0,0], then the resulting angle
+  must be PI/2.0. To get correct result, note that the common edge is AC, and
+  one must call the method as angle_btw_triangles(B,A,C,D). With this call we
+  have ensured that N1 = AB x AC (which coincides with [0,0,1]) and N2 = AC x AD
+  (which coincides with [1,0,0]). Alternatively, if the orientations of the two
+  triangles were the oppisite, the correct call would be
+  angle_btw_triangles(B,C,A,D) so that N1 = CB x CA and N2 = CA x CD.
+  */
+
+  const Utils::Vector3d a{1, 1, 1}, b{2, 1, 1}, c{1, 2, 1}, d{1, 1, 2};
+  using Utils::angle_btw_triangles;
+  BOOST_CHECK_SMALL(std::abs(angle_btw_triangles(b, a, c, d) - M_PI / 2.0),
+                    epsilon);
+  BOOST_CHECK_SMALL(std::abs(angle_btw_triangles(b, c, a, d) - 3 * M_PI / 2.0),
+                    epsilon);
+}