Oif corrections

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/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);
+};

testsuite/python/oif_volume_conservation.py

@@ -15,6 +15,7 @@
 # You should have received a copy of the GNU General Public License
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 import espressomd
+import numpy as np
 import unittest as ut
 import unittest_decorators as utx
 from tests_common import abspath
@@ -33,7 +34,6 @@ def test(self):
         self.assertEqual(system.max_oif_objects, 0)
         system.time_step = 0.4
         system.cell_system.skin = 0.5
-        system.thermostat.set_langevin(kT=0, gamma=0.7, seed=42)
 
         # creating the template for OIF object
         cell_type = oif.OifCellType(
@@ -59,10 +59,26 @@ def test(self):
         diameter_stretched = cell0.diameter()
         print("stretched diameter = " + str(diameter_stretched))
 
+        # Apply non-isotropic deformation
+        system.part[:].pos = system.part[:].pos * np.array((0.96, 1.05, 1.02))
+
+        # Test that restoring forces net to zero and don't produce a torque
+        system.integrator.run(1)
+        np.testing.assert_allclose(
+            np.sum(
+                system.part[:].f, axis=0), [
+                0., 0., 0.], atol=1E-12)
+
+        total_torque = np.zeros(3)
+        for p in system.part:
+            total_torque += np.cross(p.pos, p.f)
+        np.testing.assert_allclose(total_torque, [0., 0., 0.], atol=2E-12)
+
         # main integration loop
+        system.thermostat.set_langevin(kT=0, gamma=0.7, seed=42)
         # OIF object is let to relax into relaxed shape of the sphere
-        for _ in range(3):
-            system.integrator.run(steps=90)
+        for _ in range(2):
+            system.integrator.run(steps=240)
             diameter_final = cell0.diameter()
             print("final diameter = " + str(diameter_final))
             self.assertAlmostEqual(