Fix MA boundary term and finish demo.

div(tau) contracts with 2nd index of tau, so in the boundary term
we need tau.n, not n.tau

Now gives good results in sinatan3 demo.

demos/demo_references.bib

@@ -8,3 +8,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_ampere_3d.py

@@ -0,0 +1,94 @@
+# 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
+# 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.)
+
+
+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:20181),
+# 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
+#
+# 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()

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`
@@ -200,12 +205,8 @@ def l2_projector(self):
             bcs.append(firedrake.EquationBC(a_bc == L_bc, self._grad_phi, i))
 
             # Allow tangential movement, but only up until the end of boundary segments
-            def perp(v, n):
-                """Return component of v perpendicular to n (assumed normalized)"""
-                return v - ufl.dot(v, n) * n
-
-            a_bc = ufl.dot(perp(v_cts, n), perp(u_cts, n)) * self.ds
-            L_bc = ufl.dot(perp(v_cts, n), perp(ufl.grad(self.phi_old), n)) * 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
@@ -328,12 +329,10 @@ def equidistributor(self):
         n = ufl.FacetNormal(self.mesh)
         sigma = firedrake.TrialFunction(self.P1_ten)
         tau = firedrake.TestFunction(self.P1_ten)
-        I = ufl.Identity(self.dim)
         a = ufl.inner(tau, sigma) * self.dx
         L = (
             -ufl.dot(ufl.div(tau), ufl.grad(self.phi)) * self.dx
-            +ufl.dot(n, ufl.dot(tau, ufl.dot(I - ufl.outer(n, n), ufl.grad(self.phi))))
-            #+ (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)
@@ -459,8 +458,7 @@ def equidistributor(self):
         F = (
             ufl.inner(tau, sigma) * self.dx
             + ufl.dot(ufl.div(tau), ufl.grad(phi)) * self.dx
-            - ufl.dot(n, ufl.dot(tau, ufl.dot(I - ufl.outer(n, n), ufl.grad(phi)))) * self.ds
-            #- (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)