Fix OIF force summary function

src/python/object_in_fluid/oif_classes.py

@@ -1147,7 +1147,7 @@ def elastic_forces(
         if (el_forces[0] == 1) or (el_forces[5] == 1) or (
                 f_metric[0] == 1) or (f_metric[5] == 1):
             # initialize list
-            stretching_forces_list = np.zeros((self.mesh.points, 3))
+            stretching_forces_list = np.zeros((len(self.mesh.points), 3))
             # calculation uses edges, but results are stored for nodes
             for e in self.mesh.edges:
                 a_current_pos = e.A.get_pos()
@@ -1171,7 +1171,7 @@ def elastic_forces(
         if (el_forces[1] == 1) or (el_forces[5] == 1) or (
                 f_metric[1] == 1) or (f_metric[5] == 1):
             # initialize list
-            bending_forces_list = np.zeros((self.mesh.points, 3))
+            bending_forces_list = np.zeros((len(self.mesh.points), 3))
             # calculation uses bending incidences, but results are stored for
             # nodes
             for angle in self.mesh.angles:
@@ -1209,7 +1209,7 @@ def elastic_forces(
         if (el_forces[2] == 1) or (el_forces[5] == 1) or (
                 f_metric[2] == 1) or (f_metric[5] == 1):
             # initialize list
-            local_area_forces_list = np.zeros((self.mesh.points, 3))
+            local_area_forces_list = np.zeros((len(self.mesh.points), 3))
             # calculation uses triangles, but results are stored for nodes
             for t in self.mesh.triangles:
                 a_current_pos = t.A.get_pos()
@@ -1243,9 +1243,7 @@ def elastic_forces(
         if (el_forces[3] == 1) or (el_forces[5] == 1) or (
                 f_metric[3] == 1) or (f_metric[5] == 1):
             # initialize list
-            global_area_forces_list = []
-            for _ in self.mesh.points:
-                global_area_forces_list.append([0.0, 0.0, 0.0])
+            global_area_forces_list = np.zeros((len(self.mesh.points), 3))
             # calculation uses triangles, but results are stored for nodes
             for t in self.mesh.triangles:
                 a_current_pos = t.A.get_pos()
@@ -1275,7 +1273,7 @@ def elastic_forces(
         if (el_forces[4] == 1) or (el_forces[5] == 1) or (
                 f_metric[4] == 1) or (f_metric[5] == 1):
             # initialize list
-            volume_forces_list = np.zeros((self.mesh.points, 3))
+            volume_forces_list = np.zeros((len(self.mesh.points), 3))
             # calculation uses triangles, but results are stored for nodes
             for t in self.mesh.triangles:
                 a_current_pos = t.A.get_pos()
@@ -1344,7 +1342,7 @@ def elastic_forces(
 
         # output vtk (folded)
         if vtk_file is not None:
-            if el_forces == (0, 0, 0, 0, 0, 0):
+            if sum(el_forces) == 0:
                 raise Exception("OifCell: elastic_forces: The option elastic_forces was not used. "
                                 "Nothing to output to vtk file.")
             self.output_vtk_pos_folded(vtk_file)
@@ -1405,8 +1403,7 @@ def elastic_forces(
                     file_name=raw_data_file, data=elastic_forces_list)
 
         # return f_metric
-        if f_metric[0] + f_metric[1] + f_metric[2] + \
-                f_metric[3] + f_metric[4] + f_metric[5] > 0:
+        if sum(f_metric) > 0:
             results = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0]
             if f_metric[0] == 1:
                 results[0] = ks_f_metric

testsuite/python/oif_volume_conservation.py

@@ -27,40 +27,37 @@
 
 
 @utx.skipIfMissingFeatures("MASS")
-class OifVolumeConservation(ut.TestCase):
-
-    """Loads a soft elastic sphere via object_in_fluid, stretches it and checks
-       restoration of original volume due to elastic forces."""
+class Test(ut.TestCase):
 
     system = espressomd.System(box_l=(50., 50., 50.))
     system.time_step = 0.4
     system.cell_system.skin = 0.5
 
-    def check_relaxation(self, **kwargs):
-        self.system.part.clear()
+    def tearDown(self):
         self.system.thermostat.turn_off()
-        half_box_l = np.copy(self.system.box_l) / 2.
+        self.system.part.clear()
 
-        # creating the template for OIF object
-        cell_type = oif.OifCellType(
+    def make_cell_type(self, **kwargs):
+        return oif.OifCellType(
             nodes_file=str(tests_common.data_path("sphere393nodes.dat")),
             triangles_file=str(
                 tests_common.data_path("sphere393triangles.dat")),
             system=self.system, **kwargs, kb=1.0, kal=1.0, kag=0.1, kv=0.1,
             check_orientation=False, resize=(3.0, 3.0, 3.0))
 
-        # creating the OIF object
-        cell0 = oif.OifCell(
-            cell_type=cell_type, particle_type=0, origin=half_box_l)
+    def check_relaxation(self, **kwargs):
+        half_box_l = np.copy(self.system.box_l) / 2.
+        cell = oif.OifCell(cell_type=self.make_cell_type(**kwargs),
+                           particle_type=0, origin=half_box_l)
         partcls = self.system.part.all()
 
-        diameter_init = cell0.diameter()
+        diameter_init = cell.diameter()
         self.assertAlmostEqual(diameter_init, 24., delta=1e-3)
 
         # OIF object is being stretched by factor 1.5
         partcls.pos = (partcls.pos - half_box_l) * 1.5 + half_box_l
 
-        diameter_stretched = cell0.diameter()
+        diameter_stretched = cell.diameter()
         self.assertAlmostEqual(diameter_stretched, 36., delta=1e-3)
 
         # Apply non-isotropic deformation
@@ -85,7 +82,7 @@ def check_relaxation(self, **kwargs):
         for _ in range(20):
             self.system.integrator.run(steps=20)
             xdata.append(self.system.time)
-            ydata.append(cell0.diameter())
+            ydata.append(cell.diameter())
         # check exponential decay
         (prefactor, lam, _, diameter_final), _ = scipy.optimize.curve_fit(
             lambda x, a, b, c, d: a * np.exp(-b * x + c) + d, xdata, ydata,
@@ -94,16 +91,42 @@ def check_relaxation(self, **kwargs):
         self.assertGreater(prefactor, 0.)
         self.assertAlmostEqual(diameter_final, diameter_init, delta=0.005)
         self.assertAlmostEqual(lam, 0.0325, delta=0.0001)
+        self.system.thermostat.turn_off()
+        self.system.part.clear()
 
-    def test(self):
+    def test_00_volume_relaxation(self):
+        """
+        Load a soft elastic sphere via ``object_in_fluid``, stretch it and
+        check restoration of original volume due to elastic forces.
+        """
         self.assertEqual(self.system.max_oif_objects, 0)
-        with self.subTest(msg='linear stretching'):
+        with self.subTest(msg="linear stretching"):
             self.check_relaxation(kslin=1.)
         self.assertEqual(self.system.max_oif_objects, 1)
-        with self.subTest(msg='neo-Hookean stretching'):
+        with self.subTest(msg="neo-Hookean stretching"):
             self.check_relaxation(ks=1.)
         self.assertEqual(self.system.max_oif_objects, 2)
 
+    def test_01_elastic_forces(self):
+        half_box_l = np.copy(self.system.box_l) / 2.
+        cell = oif.OifCell(cell_type=self.make_cell_type(),
+                           particle_type=0, origin=half_box_l)
+        # stretch cell
+        partcls = self.system.part.all()
+        partcls.pos = (partcls.pos - half_box_l) * 1.5 + half_box_l
+        # reduce number of triangles to speed up calculation
+        cell.mesh.triangles = cell.mesh.triangles[:20]
+        # smoke test
+        results = cell.elastic_forces(el_forces=6 * [0], f_metric=6 * [1])
+        ref = [0., 1.36730815e-12, 22.4985704, 6838.5749, 7.3767594, 6816.6342]
+        np.testing.assert_allclose(results, ref, atol=1e-10, rtol=1e-7)
+        self.assertEqual(cell.elastic_forces(), 0)
+        # check exceptions
+        with self.assertRaises(Exception):
+            cell.elastic_forces(el_forces=6 * [0], vtk_file="test")
+        with self.assertRaises(Exception):
+            cell.elastic_forces(el_forces=6 * [0], raw_data_file="test")
+
 
 if __name__ == "__main__":
     ut.main()