Use mesh tangling checker (#103)

Closes #101.

This PR hooks up the `MeshTanglingChecker` to the mover classes and
turns it on by default in the 2D case. This will stop subtly tangled
meshes from eluding detection.

movement/laplacian_smoothing.py

@@ -99,6 +99,8 @@ def move(self, time, update_boundary_velocity=None, boundary_conditions=None):
         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"

movement/monge_ampere.py

@@ -165,6 +165,9 @@ def _update_coordinates(self):
         self.x.assign(self.xi + self.grad_phi)
         self.mesh.coordinates.assign(self.x)
 
+        if hasattr(self, "tangling_checker"):
+            self.tangling_checker.check()
+
     @property
     @PETSc.Log.EventDecorator()
     def l2_projector(self):

movement/mover.py

@@ -8,6 +8,8 @@
 from firedrake.cython.dmcommon import create_section
 from firedrake.petsc import PETSc
 
+from movement.tangling import MeshTanglingChecker
+
 __all__ = ["PrimeMover"]
 
 
@@ -17,7 +19,12 @@ class PrimeMover:
     """
 
     def __init__(
-        self, mesh, monitor_function=None, raise_convergence_errors=True, **kwargs
+        self,
+        mesh,
+        monitor_function=None,
+        raise_convergence_errors=True,
+        tangling_check=None,
+        **kwargs,
     ):
         r"""
         :arg mesh: the physical mesh
@@ -27,7 +34,10 @@ def __init__(
         :kwarg raise_convergence_errors: convergence error handling behaviour: if `True`
             then :class:`~.ConvergenceError`\s are raised, else warnings are raised and
             the program is allowed to continue
-        :kwarg raise_convergence_errors: :class:`bool`
+        :type raise_convergence_errors: :class:`bool`
+        :kwarg tangling_check: check whether the mesh has tangled elements (by default
+            on in the 2D case and off otherwise)
+        :type tangling_check: :class:`bool`
         """
         self.mesh = firedrake.Mesh(mesh.coordinates.copy(deepcopy=True))
         self.monitor_function = monitor_function
@@ -52,6 +62,14 @@ def __init__(
         self._create_function_spaces()
         self._create_functions()
 
+        # Utilities
+        if tangling_check is None:
+            tangling_check = self.dim == 2
+        if tangling_check:
+            self.tangling_checker = MeshTanglingChecker(
+                self.mesh, raise_error=raise_convergence_errors
+            )
+
     @abc.abstractmethod
     def _create_function_spaces(self):
         self.coord_space = self.mesh.coordinates.function_space()

movement/spring.py

@@ -278,6 +278,9 @@ def move(self, time, update_boundary_displacement=None, boundary_conditions=None
             f"   Displacement {np.linalg.norm(self.displacement):.2f} m"
         )
 
+        if hasattr(self, "tangling_checker"):
+            self.tangling_checker.check()
+
 
 class SpringMover_Torsional(SpringMover_Lineal):
     """

movement/tangling.py

@@ -2,6 +2,7 @@
 
 import firedrake
 import ufl
+from firedrake.__future__ import interpolate
 
 __all__ = ["MeshTanglingChecker"]
 
@@ -28,12 +29,11 @@ def __init__(self, mesh, raise_error=True):
         # Store initial signs of Jacobian determinant
         P0 = firedrake.FunctionSpace(mesh, "DG", 0)
         detJ = ufl.JacobianDeterminant(mesh)
-        s = firedrake.Function(P0)
-        s.interpolate(ufl.sign(detJ))
+        s = firedrake.assemble(interpolate(ufl.sign(detJ), P0))
 
         # Get scaled Jacobian expression
         P0_ten = firedrake.TensorFunctionSpace(mesh, "DG", 0)
-        J = firedrake.interpolate(ufl.Jacobian(mesh), P0_ten)
+        J = firedrake.assemble(interpolate(ufl.Jacobian(mesh), P0_ten))
         edge1 = ufl.as_vector([J[0, 0], J[1, 0]])
         edge2 = ufl.as_vector([J[0, 1], J[1, 1]])
         edge3 = edge1 - edge2

test/test_monge_ampere.py

@@ -35,6 +35,13 @@ def test_no_monitor_valueerror(self):
             MongeAmpereMover(self.dummy_mesh, None)
         self.assertEqual(str(cm.exception), "Please supply a monitor function.")
 
+    def test_tangling_valueerror(self):
+        mover = MongeAmpereMover(self.mesh(2, n=3), ring_monitor)
+        mover.xi.dat.data[3] += 0.2
+        with self.assertRaises(ValueError) as cm:
+            mover.move()
+        self.assertEqual(str(cm.exception), "Mesh has 1 tangled element.")
+
     @parameterized.expand(
         [(2, "relaxation"), (2, "quasi_newton"), (3, "relaxation"), (3, "quasi_newton")]
     )

test/test_spring.py

@@ -52,6 +52,13 @@ def test_convergence_error(self):
             mover.move(0)
         self.assertEqual(str(cm.exception), "Solver failed to converge.")
 
+    def test_tangling_valueerror(self):
+        mover = SpringMover(UnitSquareMesh(3, 3), 1.0)
+        mover.mesh.coordinates.dat.data[3] += 0.2
+        with self.assertRaises(ValueError) as cm:
+            mover.move(0)
+        self.assertEqual(str(cm.exception), "Mesh has 1 tangled element.")
+
 
 class TestStiffness(unittest.TestCase):
     """
@@ -191,7 +198,9 @@ def test_force_right_free(self):
     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)
+        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)