Unit tests for Laplacian smoothing (#111)

Closes #19.
Closes #78.

This will finish covering Movement with unit tests once #95 and #109 are
merged.

movement/laplacian_smoothing.py

@@ -4,7 +4,7 @@
 import ufl
 from firedrake.petsc import PETSc
 
-import movement.solver_parameters as solver_parameters
+import movement.solver_parameters as sp
 from movement.mover import PrimeMover
 
 __all__ = ["LaplacianSmoother"]
@@ -50,7 +50,15 @@ def _create_functions(self):
         self.rhs = firedrake.Function(self.coord_space, name="Zero RHS")
 
     @PETSc.Log.EventDecorator()
-    def _setup_solver(self, boundary_conditions):
+    def _solve(self, boundary_conditions, solver_parameters=None):
+        """
+        Solve the Laplace system.
+
+        :kwarg boundary_conditions: Dirichlet boundary conditions to be enforced
+        :type boundary_conditions: :class:`~.DirichletBC` or :class:`list` thereof
+        :kwarg solver_parameters: solver parameters for solving the Laplace equation
+        :type solver_parameters: :class:`dict`
+        """
         if not hasattr(self, "_solver"):
             test = firedrake.TestFunction(self.coord_space)
             trial = firedrake.TrialFunction(self.coord_space)
@@ -62,49 +70,56 @@ def _setup_solver(self, boundary_conditions):
             )
             self._solver = firedrake.LinearVariationalSolver(
                 problem,
-                solver_parameters=solver_parameters.cg_ilu,
+                solver_parameters=solver_parameters or sp.cg_ilu,
             )
         self._solver.solve()
 
     @PETSc.Log.EventDecorator()
-    def move(self, time, update_boundary_velocity=None, boundary_conditions=None):
+    def move(
+        self,
+        time,
+        update_boundary_velocity=None,
+        boundary_conditions=None,
+        solver_parameters=None,
+    ):
         """
         Assemble and solve the Laplacian system and update the coordinates.
 
         :arg time: the current time
         :type time: :class:`float`
-        :kwarg update_boundary_velocity: function that updates the boundary conditions at
-            the current time
+        :kwarg update_boundary_velocity: function that updates the boundary conditions
+            at the current time
         :type update_boundary_velocity: :class:`~.Callable` with a single argument of
             :class:`float` type
         :kwarg boundary_conditions: Dirichlet boundary conditions to be enforced
         :type boundary_conditions: :class:`~.DirichletBC` or :class:`list` thereof
+        :kwarg solver_parameters: solver parameters for solving the Laplace equation
+        :type solver_parameters: :class:`dict`
         """
         if update_boundary_velocity is not None:
             update_boundary_velocity(time)
         if not boundary_conditions:
             boundary_conditions = firedrake.DirichletBC(
                 self.coord_space, 0, "on_boundary"
             )
-        self._setup_solver(boundary_conditions)
 
         # Solve on computational mesh
         self.mesh.coordinates.assign(self.xi)
         try:
-            self._solver.solve()
+            self._solve(boundary_conditions, solver_parameters=solver_parameters)
         except fexc.ConvergenceError as conv_err:
             self._convergence_error(exception=conv_err)
 
         # Update mesh coordinates
         self.displacement[:] = self.v.dat.data_with_halos * self.dt
         self.x.dat.data_with_halos[:] += self.displacement
         self.mesh.coordinates.assign(self.x)
-        if hasattr(self, "tangling_checker"):
-            self.tangling_checker.check()
         self.volume.interpolate(ufl.CellVolume(self.mesh))
         PETSc.Sys.Print(
             f"{time:.2f} s"
             f"   Volume ratio {self.volume_ratio:5.2f}"
             f"   Variation (σ/μ) {self.coefficient_of_variation:8.2e}"
             f"   Displacement {np.linalg.norm(self.displacement):.2f} m"
         )
+        if hasattr(self, "tangling_checker"):
+            self.tangling_checker.check()

test/test_forced_movement.py

@@ -0,0 +1,151 @@
+import unittest
+
+import numpy as np
+from firedrake import *
+
+from movement import LaplacianSmoother, SpringMover
+
+
+class BaseClasses:
+    """
+    Base classes for testing mesh movement under forcings.
+    """
+
+    class TestBoundaryForcing(unittest.TestCase):
+        """
+        Unit tests for mesh movement under boundary forcings.
+        """
+
+        def move(
+            self,
+            mesh,
+            fixed_boundary_tags=None,
+            moving_boundary_tags=None,
+            vector=[1, 0],
+            **kwargs,
+        ):
+            mover = self.mover(mesh)
+            bcs = []
+            if fixed_boundary_tags:
+                bcs.append(DirichletBC(mover.coord_space, 0, fixed_boundary_tags))
+            if moving_boundary_tags:
+                f = Function(mover.coord_space).interpolate(as_vector(vector))
+                bcs.append(DirichletBC(mover.coord_space, f, moving_boundary_tags))
+            mover.move(0.0, boundary_conditions=bcs, **kwargs)
+            return mover.mesh
+
+        @staticmethod
+        def shifted_mesh(nx, ny, shift_x=0, shift_y=0):
+            mesh = UnitSquareMesh(nx, ny)
+            x, y = SpatialCoordinate(mesh)
+            shifted_coords = Function(mesh.coordinates)
+            shifted_coords.interpolate(as_vector([x + shift_x, y + shift_y]))
+            return Mesh(shifted_coords)
+
+        @staticmethod
+        def stretched_mesh(nx, ny, stretch_x=1, stretch_y=1):
+            return RectangleMesh(nx, ny, stretch_x, stretch_y)
+
+        def test_fixed_triangle(self):
+            mesh = UnitTriangleMesh()
+            coords = mesh.coordinates.dat.data
+            self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data))
+
+        def test_fixed_square(self):
+            mesh = UnitSquareMesh(1, 1)
+            coords = mesh.coordinates.dat.data
+            self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data))
+
+        def test_force_right_free(self):
+            mesh = UnitSquareMesh(10, 10)
+            coord_array = mesh.coordinates.dat.data
+            new_coord_array = self.move(
+                mesh, moving_boundary_tags=1
+            ).coordinates.dat.data
+            self.assertFalse(np.allclose(coord_array, new_coord_array))
+            shifted_coord_array = self.shifted_mesh(
+                10, 10, shift_x=1
+            ).coordinates.dat.data
+            self.assertTrue(np.allclose(shifted_coord_array, new_coord_array))
+
+        def test_force_right_fixed(self):
+            mesh = UnitSquareMesh(10, 10)
+            coord_array = mesh.coordinates.dat.data
+            new_mesh = self.move(
+                mesh, moving_boundary_tags=1, fixed_boundary_tags=2, vector=[0.1, 0]
+            )
+            new_coord_array = new_mesh.coordinates.dat.data
+            self.assertFalse(np.allclose(coord_array, new_coord_array))
+            # # TODO: Implement no-slip BCs for segments 3 and 4 (#99)
+            # stretched_mesh = self.stretched_mesh(10, 10, stretch_x=2)
+            # stretched_coord_array = stretched_mesh.coordinates.dat.data
+            # self.assertTrue(np.allclose(stretched_coord_array, new_coord_array))
+
+        def test_force_right_left_free(self):
+            mesh = UnitSquareMesh(10, 10)
+            coord_array = mesh.coordinates.dat.data
+            mesh = self.move(mesh, moving_boundary_tags=1, vector=[1, 0])
+            self.assertFalse(np.allclose(coord_array, mesh.coordinates.dat.data))
+            mesh = self.move(mesh, moving_boundary_tags=1, vector=[-1, 0])
+            self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))
+
+
+class TestSpringMover(BaseClasses.TestBoundaryForcing):
+    """
+    Unit tests for the lineal spring method under boundary forcings.
+    """
+
+    @staticmethod
+    def mover(mesh):
+        return SpringMover(mesh, 1.0)
+
+    def test_update_boundary_displacement(self):
+        mesh = UnitSquareMesh(10, 10)
+        coord_array = mesh.coordinates.dat.data
+        mover = self.mover(mesh)
+        forcing = Function(mover.coord_space)
+        bc = DirichletBC(mover.coord_space, forcing, 1)
+
+        def update_bc(time):
+            forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)]))
+
+        mover.move(0.0, boundary_conditions=bc, update_boundary_displacement=update_bc)
+        self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data))
+        mover.move(1.0, boundary_conditions=bc, update_boundary_displacement=update_bc)
+        self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))
+
+
+class TestLaplacianSmoother(BaseClasses.TestBoundaryForcing):
+    """
+    Unit tests for Laplacian smoothing under boundary forcings.
+    """
+
+    @staticmethod
+    def mover(mesh):
+        return LaplacianSmoother(mesh, 1.0)
+
+    def test_convergence_error(self):
+        with self.assertRaises(ConvergenceError) as cm:
+            self.move(
+                UnitSquareMesh(10, 10),
+                moving_boundary_tags=1,
+                fixed_boundary_tags=2,
+                vector=[1.0, 0],
+                solver_parameters={"ksp_type": "cg", "ksp_max_it": 0},
+            )
+        self.assertEqual(str(cm.exception), "Solver failed to converge.")
+
+    def test_update_boundary_velocity(self):
+        mesh = UnitSquareMesh(10, 10)
+        coord_array = mesh.coordinates.dat.data
+        mover = LaplacianSmoother(mesh, 1.0)
+        forcing = Function(mover.coord_space)
+        bc = DirichletBC(mover.coord_space, forcing, 1)
+
+        def update_bc(time):
+            forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)]))
+
+        mover.move(0.0, boundary_conditions=bc, update_boundary_velocity=update_bc)
+        self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data))
+        mover.move(1.0, boundary_conditions=bc, update_boundary_velocity=update_bc)
+        self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))

test/test_laplacian_smoothing.py

@@ -0,0 +1,20 @@
+import unittest
+from unittest.mock import MagicMock
+
+from firedrake import *
+
+from movement import LaplacianSmoother
+
+
+class TestExceptions(unittest.TestCase):
+    """
+    Unit tests for exceptions raised during Laplacian smoothing.
+    """
+
+    def test_tangling_valueerror(self):
+        mover = LaplacianSmoother(UnitSquareMesh(3, 3), 1.0)
+        mover.x.dat.data[3] += 0.2
+        mover._solve = MagicMock()
+        with self.assertRaises(ValueError) as cm:
+            mover.move(0)
+        self.assertEqual(str(cm.exception), "Mesh has 1 tangled element.")

test/test_spring.py

@@ -146,87 +146,3 @@ def test_angles(self):
         """
         angles = self.rad2deg(self.mover.angles.dat.data)
         self.assertTrue(np.allclose(angles, [0, 135, 90]))
-
-
-class TestMovement(unittest.TestCase):
-    """
-    Unit tests for movement under spring-based methods.
-    """
-
-    @staticmethod
-    def move(mesh, fixed_boundary_tags=None, moving_boundary_tags=None, vector=[1, 0]):
-        mover = SpringMover(mesh, 1.0)
-        bcs = []
-        if fixed_boundary_tags:
-            bcs.append(DirichletBC(mover.coord_space, 0, fixed_boundary_tags))
-        if moving_boundary_tags:
-            f = Function(mover.coord_space).interpolate(as_vector(vector))
-            bcs.append(DirichletBC(mover.coord_space, f, moving_boundary_tags))
-        mover.move(0.0, boundary_conditions=bcs)
-        return mover.mesh
-
-    @staticmethod
-    def shifted_mesh(nx, ny, shift_x=0, shift_y=0):
-        mesh = UnitSquareMesh(nx, ny)
-        x, y = SpatialCoordinate(mesh)
-        shifted_coords = Function(mesh.coordinates)
-        shifted_coords.interpolate(as_vector([x + shift_x, y + shift_y]))
-        return Mesh(shifted_coords)
-
-    @staticmethod
-    def stretched_mesh(nx, ny, stretch_x=1, stretch_y=1):
-        return RectangleMesh(nx, ny, stretch_x, stretch_y)
-
-    def test_fixed_triangle(self):
-        mesh = UnitTriangleMesh()
-        coords = mesh.coordinates.dat.data
-        self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data))
-
-    def test_fixed_square(self):
-        mesh = UnitSquareMesh(1, 1)
-        coords = mesh.coordinates.dat.data
-        self.assertTrue(np.allclose(coords, self.move(mesh).coordinates.dat.data))
-
-    def test_force_right_free(self):
-        mesh = UnitSquareMesh(10, 10)
-        coord_array = mesh.coordinates.dat.data
-        new_coord_array = self.move(mesh, moving_boundary_tags=1).coordinates.dat.data
-        self.assertFalse(np.allclose(coord_array, new_coord_array))
-        shifted_coord_array = self.shifted_mesh(10, 10, shift_x=1).coordinates.dat.data
-        self.assertTrue(np.allclose(shifted_coord_array, new_coord_array))
-
-    def test_force_right_fixed(self):
-        mesh = UnitSquareMesh(10, 10)
-        coord_array = mesh.coordinates.dat.data
-        new_mesh = self.move(
-            mesh, moving_boundary_tags=1, fixed_boundary_tags=2, vector=[0.1, 0]
-        )
-        new_coord_array = new_mesh.coordinates.dat.data
-        self.assertFalse(np.allclose(coord_array, new_coord_array))
-        # # TODO: Implement no-slip BCs for segments 3 and 4 (#99)
-        # stretched_mesh = self.stretched_mesh(10, 10, stretch_x=2)
-        # stretched_coord_array = stretched_mesh.coordinates.dat.data
-        # self.assertTrue(np.allclose(stretched_coord_array, new_coord_array))
-
-    def test_force_right_left_free(self):
-        mesh = UnitSquareMesh(10, 10)
-        coord_array = mesh.coordinates.dat.data
-        mesh = self.move(mesh, moving_boundary_tags=1, vector=[1, 0])
-        self.assertFalse(np.allclose(coord_array, mesh.coordinates.dat.data))
-        mesh = self.move(mesh, moving_boundary_tags=1, vector=[-1, 0])
-        self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))
-
-    def test_update_boundary_displacement(self):
-        mesh = UnitSquareMesh(10, 10)
-        coord_array = mesh.coordinates.dat.data
-        mover = SpringMover(mesh, 1.0)
-        forcing = Function(mover.coord_space)
-        bc = DirichletBC(mover.coord_space, forcing, 1)
-
-        def update_bc(time):
-            forcing.interpolate(as_vector([cos(pi * time / 2), sin(pi * time / 2)]))
-
-        mover.move(0.0, boundary_conditions=bc, update_boundary_displacement=update_bc)
-        self.assertFalse(np.allclose(coord_array, mover.mesh.coordinates.dat.data))
-        mover.move(1.0, boundary_conditions=bc, update_boundary_displacement=update_bc)
-        self.assertTrue(np.allclose(coord_array, mesh.coordinates.dat.data))