From fb817ad82ee2cadebdb3b825cdb0c604c8a94551 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 20 Oct 2023 15:04:30 +0200 Subject: [PATCH 1/2] OIF: Repair OifCell.elastic_forces() feature Fix regression introduced in 56e77ed7547faea8fdc293278817e0addce769dd. --- src/python/object_in_fluid/oif_classes.py | 17 +++++++---------- 1 file changed, 7 insertions(+), 10 deletions(-) diff --git a/src/python/object_in_fluid/oif_classes.py b/src/python/object_in_fluid/oif_classes.py index f471f18d45f..7d5fa4e956b 100644 --- a/src/python/object_in_fluid/oif_classes.py +++ b/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 From 4ae85fdf477df2165a2b4c2a5b981385df6abee4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jean-No=C3=ABl=20Grad?= Date: Fri, 20 Oct 2023 16:00:42 +0200 Subject: [PATCH 2/2] tests: Add OIF smoke test --- testsuite/python/oif_volume_conservation.py | 59 ++++++++++++++------- 1 file changed, 41 insertions(+), 18 deletions(-) diff --git a/testsuite/python/oif_volume_conservation.py b/testsuite/python/oif_volume_conservation.py index 054787f24f0..ca1b0673e7a 100644 --- a/testsuite/python/oif_volume_conservation.py +++ b/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()