From c6b51c100c8840fa5a8aca18ed44ca7504369a94 Mon Sep 17 00:00:00 2001
From: Joe Wallwork <22053413+jwallwork23@users.noreply.github.com>
Date: Thu, 4 Apr 2024 17:14:51 +0100
Subject: [PATCH] Add Monge-Ampere zoom plot and tangling exercise (#65)

Closes #50. Partly addresses #19.

This PR adds a zoom plot to the first Monge-Ampere demo to show that the
elements aren't tangled, as well as an exercise to check this with
`MeshTanglingChecker`.

The PR also applies some minor refactoring to `MeshTanglingChecker`, as
well as writing unit tests for it.
---
 demos/monge_ampere1.py | 21 +++++++++++++++++++++
 movement/tangling.py   | 41 +++++++++++++++++++----------------------
 test/test_tangling.py  | 37 ++++++++++++++++++++++++++++---------
 3 files changed, 68 insertions(+), 31 deletions(-)

diff --git a/demos/monge_ampere1.py b/demos/monge_ampere1.py
index d3fd2e3..dcae8cf 100644
--- a/demos/monge_ampere1.py
+++ b/demos/monge_ampere1.py
@@ -103,6 +103,27 @@ def ring_monitor(mesh):
 #    :figwidth: 60%
 #    :align: center
 #
+# Whilst it looks like the mesh might have tangled elements, closer inspection shows
+# that this is not the case.
+
+fig, axes = plt.subplots()
+triplot(mover.mesh, axes=axes)
+axes.set_xlim([0.15, 0.3])
+axes.set_ylim([0.15, 0.3])
+axes.set_aspect(1)
+plt.savefig("monge_ampere1-adapted_mesh_zoom.jpg")
+
+# .. figure:: monge_ampere1-adapted_mesh_zoom.jpg
+#    :figwidth: 60%
+#    :align: center
+#
+# .. rubric:: Exercise
+#
+# To further convince yourself that there are no tangled elements, go back to the start
+# of the demo and set up a :class:`movement.tangling.MeshTanglingChecker` object using
+# the initial mesh. Use it to check for tangling after the mesh movement has been
+# applied.
+#
 # This tutorial can be dowloaded as a `Python script <monge_ampere1.py>`__.
 #
 #
diff --git a/movement/tangling.py b/movement/tangling.py
index 5c14c77..f9f9064 100644
--- a/movement/tangling.py
+++ b/movement/tangling.py
@@ -1,30 +1,29 @@
 import firedrake
 import ufl
+import warnings
 
 
 __all__ = ["MeshTanglingChecker"]
 
 
-class MeshTanglingChecker(object):
+class MeshTanglingChecker:
     """
-    A class for tracking whether a mesh has
-    tangled, i.e. whether any of its elements
+    A class for tracking whether a mesh has tangled, i.e. whether any of its elements
     have become inverted.
     """
 
-    def __init__(self, mesh, mode="warn"):
+    def __init__(self, mesh, raise_error=True):
         """
         :arg mesh: the mesh to track if tangled
-        :kwarg mode: should a warning or an error
-            be raised when tangling is encountered?
+        :type mesh: :class:`firedrake.mesh.MeshGeometry`
+        :kwarg raise_error: if ``True``, an error is raised if any element is tangled,
+            otherwise a warning is raised
+        :type raise_error: :class:`bool`
         """
-        dim = mesh.topological_dimension()
-        if dim != 2:
-            raise ValueError(f"Cannot check for tangling of {dim}D meshes (only 2D)")
+        if mesh.topological_dimension() != 2:
+            raise NotImplementedError("Tangling check only currently supported in 2D.")
         self.mesh = mesh
-        if mode not in ["warn", "error"]:
-            raise ValueError("Choose mode from 'warn' and 'error'")
-        self.mode = mode
+        self.raise_error = raise_error
 
         # Store initial signs of Jacobian determinant
         P0 = firedrake.FunctionSpace(mesh, "DG", 0)
@@ -55,17 +54,15 @@ def scaled_jacobian(self):
 
     def check(self):
         """
-        Check for tangling.
+        Check whether any element orientations have changed since the tangling checker
+        was created.
         """
-        sj = self.scaled_jacobian.dat.data
+        sj = self.scaled_jacobian.dat.data_with_halos
         num_tangled = len(sj[sj < 0])
-        if num_tangled == 0:
-            return 0
-        msg = f"Mesh has {num_tangled} tangled elements"
-        if self.mode == "warn":
-            import warnings
-
+        if num_tangled > 0:
+            plural = "s" if num_tangled > 1 else ""
+            msg = f"Mesh has {num_tangled} tangled element{plural}."
+            if self.raise_error:
+                raise ValueError(msg)
             warnings.warn(msg)
-        else:
-            raise ValueError(msg)
         return num_tangled
diff --git a/test/test_tangling.py b/test/test_tangling.py
index 59b56e8..a5fe278 100644
--- a/test/test_tangling.py
+++ b/test/test_tangling.py
@@ -1,15 +1,34 @@
 from firedrake import *
 from movement import *
+import unittest
 
 
-def test_tangling_checker():
+class TestTangling(unittest.TestCase):
     """
-    Test that :class:`MeshTanglingChecker`
-    can correctly identify tangled elements
-    in a 2D triangular mesh.
+    Unit tests for mesh tangling checking.
     """
-    mesh = UnitSquareMesh(3, 3)
-    checker = MeshTanglingChecker(mesh)
-    assert checker.check() == 0
-    mesh.coordinates.dat.data[3] += 0.2
-    assert checker.check() == 1
+
+    def test_tangling_checker_error1(self):
+        mesh = UnitSquareMesh(3, 3)
+        checker = MeshTanglingChecker(mesh, raise_error=True)
+        mesh.coordinates.dat.data[3] += 0.2
+        with self.assertRaises(ValueError) as cm:
+            checker.check()
+        msg = "Mesh has 1 tangled element."
+        self.assertEqual(str(cm.exception), msg)
+
+    def test_tangling_checker_error2(self):
+        mesh = UnitSquareMesh(3, 3)
+        checker = MeshTanglingChecker(mesh, raise_error=True)
+        mesh.coordinates.dat.data[3] += 0.5
+        with self.assertRaises(ValueError) as cm:
+            checker.check()
+        msg = "Mesh has 3 tangled elements."
+        self.assertEqual(str(cm.exception), msg)
+
+    def test_tangling_checker_warning1(self):
+        mesh = UnitSquareMesh(3, 3)
+        checker = MeshTanglingChecker(mesh, raise_error=False)
+        self.assertEqual(checker.check(), 0)
+        mesh.coordinates.dat.data[3] += 0.2
+        self.assertEqual(checker.check(), 1)