Extend MA mover to 3D and add demo

demos/demo_references.bib

@@ -19,3 +19,12 @@ @Article{McRae:2018
   pages={1121--1148},
   doi={10.1137/16M1109515}
 }
+@inproceedings{park2019,
+  title={Verification of unstructured grid adaptation components},
+  author={Park, Michael A and Balan, Aravind and Anderson, William K and Galbraith, Marshall C and Caplan, Philip and Carson, Hugh A and Michal, Todd R and Krakos, Joshua A and Kamenetskiy, Dmitry S and Loseille, Adrien and others},
+  booktitle={AIAA Scitech 2019 Forum},
+  pages={1723},
+  year={2019},
+  doi={10.2514/6.2019-1723}
+}
+

demos/monge_ampere1.py

@@ -124,6 +124,9 @@ def ring_monitor(mesh):
 # the initial mesh. Use it to check for tangling after the mesh movement has been
 # applied.
 #
+# In the `next demo <./monge_ampere_3d.py.html>`__, we will demonstrate
+# that the Monge Ampère method can also be applied in three dimensions.
+#
 # This tutorial can be dowloaded as a `Python script <monge_ampere1.py>`__.
 #
 #

demos/monge_ampere_3d.py

@@ -0,0 +1,99 @@
+# Monge-Ampère mesh movement in three dimensions
+# ==============================================
+
+# In this demo we demonstrate that the  Monge-Ampère mesh movement
+# can also be applied to 3D meshes. We employ the `sinatan3` function
+# from :cite:`park2019` to introduce an interesting pattern.
+
+import movement
+from firedrake import *
+
+
+def sinatan3(mesh):
+    x, y, z = SpatialCoordinate(mesh)
+    return 0.1 * sin(50 * x * z) + atan2(0.1, sin(5 * y) - 2 * x * z)
+
+
+# We will now try to use mesh movement to optimize the mesh such that it can
+# most accurately represent this function with limited resolution.
+# A good indicator for where resolution is required
+# is to look at the curvature of the function which can be expressed
+# in terms of the norm of the Hessian. A monitor function
+# that targets high resolution in places of high curvature then looks like
+#
+# .. math::
+#
+#    m = 1 + \alpha \frac{H(u_h):H(u_h)}{\max_{{\bf x}\in\Omega} H(u_h):H(u_h)}
+#
+# where :math:`:` indicates the inner product, i.e. :math:`\sqrt{H:H}` is the Frobenius norm
+# of :math:`H`. We have normalized such that the minimum of the monitor function is one (where
+# the error is zero), and its maximum is :math:`1 + \alpha` (where the curvature is maximal). This
+# means we can select the ratio between the largest and smallest cell volume  in the
+# moved mesh as :math:`1+\alpha`.
+#
+# As in the `previous Monge-Ampère demo <./monge_ampere1.py.html>`__, we use the
+# :class:`~.MongeAmpereMover` to perform the mesh movement based on this monitor. We need
+# to provide the monitor as a callback function that takes the mesh as its
+# input. During the iterations of the mesh movement process the monitor will then
+# be re-evaluated in the (iteratively)
+# moved mesh nodes so that, as we improve the mesh, we can also more accurately
+# express the monitor function in the desired high-resolution areas.
+# To track what happens during these iterations, we define a VTK file object
+# that we will write to in every call when the monitor gets re-evaluated.
+
+from firedrake.output import VTKFile
+
+f = VTKFile("monitor.pvd")
+alpha = Constant(10.0)
+
+
+def monitor(mesh):
+    V = FunctionSpace(mesh, "CG", 1)
+    # interpolation of the function itself, for output purposes only:
+    u = Function(V, name="sinatan3")
+    u.interpolate(sinatan3(mesh))
+
+    # NOTE: we are taking the Hessian of the _symbolic_ UFL expression
+    # returned by sinatan3(mesh), *not* of the P1 interpolated version
+    # stored in u. u is a piecewise linear approximation; Taking its gradient
+    # once would result in a piecewise constant (cell-wise) gradient, and taking
+    # the gradient of that again would simply be zero.
+    hess = grad(grad(sinatan3(mesh)))
+
+    hnorm = Function(V, name="hessian_norm")
+    hnorm.interpolate(inner(hess, hess))
+    hmax = hnorm.dat.data[:].max()
+    f.write(u, hnorm)
+    return 1.0 + alpha * hnorm / Constant(hmax)
+
+
+# Let us try this on a fairly coarse unit cube mesh. Note that
+# we use the `"relaxation"` method (see :cite:`McRae:2018`),
+# which gives faster convergence for this case.
+
+n = 20
+mesh = UnitCubeMesh(n, n, n)
+mover = movement.MongeAmpereMover(mesh, monitor, method="relaxation")
+mover.move()
+
+# The results will be written to the `monitor.pvd` file which represents
+# a series of outputs storing the function and hessian norm evaluated
+# on each of the iteratively moved meshes. This can be viewed, e.g., using
+# `ParaView <https://www.paraview.org/>`_, to produce the following
+# image:
+#
+# .. figure:: monge_ampere_3d-paraview.jpg
+#    :figwidth: 60%
+#    :align: center
+#
+# In the `next demo <./monge_ampere_helmholtz.py.html>`__, we will demonstrate
+# how to optimize the mesh for the discretisation of a PDE with the aim to
+# minimize its discretisation error.
+#
+# This tutorial can be dowloaded as a `Python script <monge_ampere_3d.py>`__.
+#
+#
+# .. rubric:: References
+#
+# .. bibliography:: demo_references.bib
+#    :filter: docname in docnames

demos/monge_ampere_helmholtz.py

@@ -96,21 +96,9 @@ def solve_helmholtz(mesh):
 #    L2-norm error on initial mesh: 0.010233816824277465
 #
 # We will now try to use mesh movement to optimize the mesh to reduce
-# this numerical error. A good indicator for where resolution is required
-# is to look at the curvature of the solution which can be expressed
-# in terms of the norm of the Hessian. A monitor function
-# that targets high resolution in places of high curvature then looks like
-#
-# .. math::
-#
-#    m = 1 + \alpha \frac{H(u_h):H(u_h)}{\max_{{\bf x}\in\Omega} H(u_h):H(u_h)}
-#
-# where :math:`:` indicates the inner product, i.e. :math:`\sqrt{H:H}` is the Frobenius norm
-# of :math:`H`. We have normalized such that the minimum of the monitor function is one (where
-# the error is zero), and its maximum is :math:`1 + \alpha` (where the curvature is maximal). This
-# means we can select the ratio between the largest area and smallest area triangle in the
-# moved mesh as :math:`1+\alpha`.
-#
+# this numerical error. We use the same monitor function as
+# in the `previous Monge-Ampère demo <./monge_ampere_3d.py.html>`__
+# based on the norm of the Hessian of the solution.
 # In the following implementation we use the exact solution :math:`u_{\text{exact}}` which we
 # have as a symbolic UFL expression, and thus we can also obtain the Hessian symbolically as
 # :code:`grad(grad(u_exact))`. To compute its maximum norm however we do interpolate it
@@ -129,6 +117,8 @@ def monitor(mesh):
     return m
 
 
+# Plot the monitor function on the original mesh
+
 fig, axes = plt.subplots()
 m = Function(u_h, name="monitor")
 m.interpolate(monitor(mesh))
@@ -140,14 +130,6 @@ def monitor(mesh):
 # .. figure:: monge_ampere_helmholtz-monitor.jpg
 #    :figwidth: 60%
 #    :align: center
-#
-# As in the `previous Monge-Ampère demo <./monge_ampere1.py.html>`__, we use the
-# MongeAmpereMover to perform the mesh movement based on this monitor. We need
-# to provide the monitor as a callback function that takes the mesh as its
-# input just as we have defined it above. During the iterations of the mesh
-# movement process the monitor will then be re-evaluated in the (iteratively)
-# moved mesh nodes so that, as we improve the mesh, we can also more accurately
-# express the monitor function in the desired high-resolution areas.
 
 mover = MongeAmpereMover(mesh, monitor, method="quasi_newton")
 mover.move()
@@ -246,3 +228,5 @@ def monitor2(mesh):
 # .. code-block:: none
 #
 #    L2-norm error on moved mesh: 0.00630874419681285
+#
+# This tutorial can be dowloaded as a `Python script <monge_ampere_helmholtz.py>`__.

movement/monge_ampere.py

@@ -15,6 +15,11 @@
 ]
 
 
+def tangential(v, n):
+    """Return component of v perpendicular to n (assumed normalized)"""
+    return v - ufl.dot(v, n) * n
+
+
 def MongeAmpereMover(mesh, monitor_function, method="relaxation", **kwargs):
     r"""
         Movement of a `mesh` is determined by a `monitor_function`
@@ -185,27 +190,23 @@ def l2_projector(self):
 
             # Check for axis-aligned boundaries
             _n = [firedrake.assemble(abs(n[j]) * self.ds(i)) for j in range(self.dim)]
-            if np.allclose(_n, 0.0):
+            iszero = [np.allclose(ni, 0.0) for ni in _n]
+            nzero = sum(iszero)
+            if nzero == self.dim:
                 raise ValueError(f"Invalid normal vector {_n}")
-            else:
-                if self.dim != 2:
-                    raise NotImplementedError  # TODO
-                if np.isclose(_n[0], 0.0):
-                    bcs.append(firedrake.DirichletBC(self.P1_vec.sub(1), 0, i))
-                    continue
-                elif np.isclose(_n[1], 0.0):
-                    bcs.append(firedrake.DirichletBC(self.P1_vec.sub(0), 0, i))
-                    continue
+            elif nzero == self.dim-1:
+                idx = iszero.index(False)
+                bcs.append(firedrake.DirichletBC(self.P1_vec.sub(idx), 0, i))
+                continue
 
             # Enforce no mesh movement normal to boundaries
             a_bc = ufl.dot(v_cts, n) * ufl.dot(u_cts, n) * self.ds
             L_bc = ufl.dot(v_cts, n) * firedrake.Constant(0.0) * self.ds
             bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, i))
 
             # Allow tangential movement, but only up until the end of boundary segments
-            s = ufl.perp(n)
-            a_bc = ufl.dot(v_cts, s) * ufl.dot(u_cts, s) * self.ds
-            L_bc = ufl.dot(v_cts, s) * ufl.dot(ufl.grad(self.phi_old), s) * self.ds
+            a_bc = ufl.dot(tangential(v_cts, n), tangential(u_cts, n)) * self.ds
+            L_bc = ufl.dot(tangential(v_cts, n), tangential(ufl.grad(self.phi_old), n)) * self.ds
             edges = set(self.mesh.exterior_facets.unique_markers)
             if len(edges) == 0:
                 bbc = None  # Periodic case
@@ -325,15 +326,13 @@ def equidistributor(self):
             return self._equidistributor
         assert hasattr(self, "phi")
         assert hasattr(self, "sigma")
-        if self.dim != 2:
-            raise NotImplementedError  # TODO
         n = ufl.FacetNormal(self.mesh)
         sigma = firedrake.TrialFunction(self.P1_ten)
         tau = firedrake.TestFunction(self.P1_ten)
         a = ufl.inner(tau, sigma) * self.dx
         L = (
             -ufl.dot(ufl.div(tau), ufl.grad(self.phi)) * self.dx
-            + (tau[0, 1] * n[1] * self.phi.dx(0) + tau[1, 0] * n[0] * self.phi.dx(1))
+            + ufl.dot(ufl.dot(tangential(ufl.grad(self.phi), n), tau), n)
             * self.ds
         )
         problem = firedrake.LinearVariationalProblem(a, L, self.sigma)
@@ -452,16 +451,14 @@ def equidistributor(self):
         if hasattr(self, "_equidistributor"):
             return self._equidistributor
         assert hasattr(self, "phisigma")
-        if self.dim != 2:
-            raise NotImplementedError  # TODO
         n = ufl.FacetNormal(self.mesh)
         I = ufl.Identity(self.dim)
         phi, sigma = firedrake.split(self.phisigma)
         psi, tau = firedrake.TestFunctions(self.V)
         F = (
             ufl.inner(tau, sigma) * self.dx
             + ufl.dot(ufl.div(tau), ufl.grad(phi)) * self.dx
-            - (tau[0, 1] * n[1] * phi.dx(0) + tau[1, 0] * n[0] * phi.dx(1)) * self.ds
+            - ufl.dot(ufl.dot(tangential(ufl.grad(phi), n), tau), n) * self.ds
             - psi * (self.monitor * ufl.det(I + sigma) - self.theta) * self.dx
         )
         phi, sigma = firedrake.TrialFunctions(self.V)

movement/mover.py

@@ -15,10 +15,6 @@ def __init__(self, mesh, monitor_function=None, **kwargs):
         self.mesh = firedrake.Mesh(mesh.coordinates.copy(deepcopy=True))
         self.monitor_function = monitor_function
         self.dim = self.mesh.topological_dimension()
-        if self.dim != 2:
-            raise NotImplementedError(
-                f"Dimension {self.dim} has not been considered yet"
-            )
         self.gdim = self.mesh.geometric_dimension()
         self.plex = self.mesh.topology_dm
         self.vertex_indices = self.plex.getDepthStratum(0)

test/monitors.py

@@ -9,6 +9,7 @@ def ring_monitor(mesh):
     alpha = Constant(10.0)  # amplitude
     beta = Constant(200.0)  # width
     gamma = Constant(0.15)  # radius
-    x, y = SpatialCoordinate(mesh)
-    r = (x - 0.5) ** 2 + (y - 0.5) ** 2
+    dim = mesh.geometric_dimension()
+    xyz = SpatialCoordinate(mesh) - as_vector([0.5]*dim)
+    r = dot(xyz, xyz)
     return Constant(1.0) + alpha / cosh(beta * (r - gamma)) ** 2

test/test_monge_ampere.py

@@ -14,13 +14,25 @@ def fix_boundary(request):
     return request.param
 
 
-def test_uniform_monitor(method, exports=False):
+@pytest.fixture(params=[2, 3])
+def dimension(request):
+    return request.param
+
+
+def uniform_mesh(dim, n):
+    if dim == 2:
+        return UnitSquareMesh(n, n)
+    else:
+        return UnitCubeMesh(n, n, n)
+
+
+def test_uniform_monitor(dimension, method, exports=False):
     """
     Test that the mesh mover converges in one
     iteration for a constant monitor function.
     """
     n = 10
-    mesh = UnitSquareMesh(n, n)
+    mesh = uniform_mesh(dimension, n)
     coords = mesh.coordinates.dat.data.copy()
 
     mover = MongeAmpereMover(mesh, const_monitor, method=method)
@@ -30,13 +42,13 @@ def test_uniform_monitor(method, exports=False):
     assert num_iterations == 0
 
 
-def test_continue(method, exports=False):
+def test_continue(dimension, method, exports=False):
     """
     Test that providing a good initial guess
     benefits the solver.
     """
     n = 20
-    mesh = UnitSquareMesh(n, n)
+    mesh = uniform_mesh(dimension, n)
     rtol = 1.0e-03
 
     # Solve the problem to a weak tolerance
@@ -52,7 +64,7 @@ def test_continue(method, exports=False):
         File("outputs/continue.pvd").write(mover.phi, mover.sigma)
 
     # Solve the problem again to a tight tolerance
-    mesh = UnitSquareMesh(n, n)
+    mesh = uniform_mesh(dimension, n)
     mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=rtol)
     num_it_naive = mover.move()
     if exports:
@@ -64,19 +76,20 @@ def test_continue(method, exports=False):
     #        for the relaxation method, which is concerning.
 
 
-def test_change_monitor(method, exports=False):
+def test_change_monitor(dimension, method, exports=False):
     """
     Test that the mover can handle changes to
     the monitor function, such as would happen
     during timestepping.
     """
     n = 20
-    mesh = UnitSquareMesh(n, n)
+    mesh = uniform_mesh(dimension, n)
+    dim = mesh.geometric_dimension()
     coords = mesh.coordinates.dat.data.copy()
     tol = 1.0e-03
 
     # Adapt to a ring monitor
-    mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=tol)
+    mover = MongeAmpereMover(mesh, ring_monitor, method=method, rtol=1e3 * tol**dim)
     mover.move()
     if exports:
         File("outputs/ring.pvd").write(mover.phi, mover.sigma)
@@ -91,13 +104,13 @@ def test_change_monitor(method, exports=False):
 
 
 @pytest.mark.slow
-def test_bcs(method, fix_boundary):
+def test_bcs(dimension, method, fix_boundary):
     """
     Test that domain boundaries are fixed by
     the Monge-Ampere movers.
     """
     n = 20
-    mesh = UnitSquareMesh(n, n)
+    mesh = uniform_mesh(dimension, n)
     one = Constant(1.0)
     bnd = assemble(one * ds(domain=mesh))
     bnodes = DirichletBC(mesh.coordinates.function_space(), 0, "on_boundary").nodes