Add basic Monge-Ampere demo

.flake8

@@ -0,0 +1,18 @@
+[flake8]
+ignore =
+    # Line too long
+    E501,
+    # Missing whitespace around arithmetic operator
+    E226,
+    # Module level import not at top of file
+    E402,
+    # Do not use variables named "I", "O", or "l"
+    E741,
+    # Unable to detect undefined names
+    F403,
+    # Name may be undefined, or defined from star imports
+    F405,
+    # Line break occurred before binary operator
+    W503
+max-line-length = 88
+select = B,C,E,F,W,T4,B9

.gitignore

@@ -1,7 +1,6 @@
-plots/
 __pycache__/
 movement.egg-info/
-*.ipynb_checkpoints/
 
+*.jpg
 *.pvd
 *.vtu

Makefile

@@ -12,7 +12,7 @@ install:
 
 lint:
 	@echo "Checking lint..."
-	@flake8 --ignore=E501,E226,E402,E731,E741,F403,F405,F999,N803,N806,W503
+	@flake8
 	@echo "PASS"
 
 test: lint

demos/demo_references.bib

@@ -0,0 +1,10 @@
+@Article{McRae:2018,
+  author={McRae, A. T. T. and Cotter, C. J. and Budd, C., J.},
+  title={Optimal-transport-based mesh adaptivity on the plane and sphere using finite elements},
+  journal={SIAM Journal on Scientific Computing},
+  year={2018},
+  volume={40},
+  number={2},
+  pages={1121--1148},
+  doi={10.1137/16M1109515}
+}

demos/monge_ampere1.py

@@ -0,0 +1,112 @@
+# Introduction to mesh movement driven by a Monge-Ampère type equation
+# ====================================================================
+
+# In this demo, we consider an example of mesh movement driven by solutions of a
+# Monge-Ampère type equation.
+#
+# The idea is to consider two domains: the *physical domain* :math:`\Omega_P` and the
+# *computational domain* :math:`\Omega_C`. Associated with these are the *physical mesh*
+# :math:`\mathcal{H}_P` and the *computational mesh* :math:`\mathcal{H}_C`. The
+# computational domain and mesh remain fixed throughout the simulation, whereas the
+# physical domain and mesh are allowed to change. In this framework, we can interpret
+# mesh movement algorithms as searches for mappings between the computational and
+# physical domains:
+#
+# .. math::
+#     \mathbf{x}:\Omega_C\rightarrow\Omega_P.
+#
+# In practice we are really searching for a discrete mapping between the computational
+# and physical meshes.
+#
+# Let :math:`\boldsymbol{\xi}` and :math:`\mathbf{x}` denote the coordinate fields in
+# the computational and physical domains. Skipping some of the details that can be
+# found in :cite:`McRae:2018`, of the possible mappings we choose one that takes the form
+#
+# .. math::
+#     \mathbf{x} = \boldsymbol{\xi}+\nabla_{\boldsymbol{\xi}}\phi,
+#
+# where :math:`\phi` is a convex potential function. Further, we choose the potential
+# such that is the solution of the Monge-Ampère type equation,
+#
+# .. math::
+#     m(\mathbf{x}) \det(\underline{\mathbf{I}}+\nabla_{\boldsymbol{\xi}}\nabla_{\boldsymbol{\xi}}\phi) = \theta,
+#
+# where :math:`m(\mathbf{x})` is a so-called *monitor function* and :math:`\theta` is a
+# strictly positive normalisation function. The monitor function is of key importance
+# because it specifies the desired *density* of the adapted mesh across the domain,
+# i.e., where resolution is focused. (Note that density is the reciprocal of area in 2D
+# or of volume in 3D.)
+#
+# We begin the example by importing from the namespaces of Firedrake and Movement.
+
+from firedrake import *
+from movement import *
+
+# To start with a simple example, consider a uniform mesh of the unit square.
+
+n = 20
+mesh = UnitSquareMesh(n, n)
+
+# We can plot the initial mesh using Matplotlib as follows.
+
+import matplotlib.pyplot as plt
+from firedrake.pyplot import triplot
+
+fig, axes = plt.subplots()
+triplot(mesh, axes=axes)
+axes.set_aspect(1)
+plt.savefig("monge_ampere1-initial_mesh.jpg")
+
+# .. figure:: monge_ampere1-initial_mesh.jpg
+#    :figwidth: 60%
+#    :align: center
+#
+# Let's choose a monitor function which focuses mesh density in a ring centred within
+# the domain, according to the formula,
+#
+# .. math::
+#     m(x,y) = 1 + \frac{\alpha}{\cosh^2\left(\beta\left(\left(x-\frac{1}{2}\right)^2+\left(y-\frac{1}{2}\right)^2-\gamma\right)\right)},
+#
+# for some values of the parameters :math:`\alpha`, :math:`\beta`, and :math:`\gamma`.
+# Unity is added at the start to ensure that the monitor function doesn't get too
+# close to zero.
+#
+# Here we can think of :math:`\alpha` as relating to the amplitude of the monitor
+# function, :math:`\beta` as relating to the width of the ring, and :math:`\gamma` as
+# the radius of the ring.
+
+
+def ring_monitor(mesh):
+    alpha = Constant(20.0)
+    beta = Constant(200.0)
+    gamma = Constant(0.15)
+    x, y = SpatialCoordinate(mesh)
+    r = (x - 0.5) ** 2 + (y - 0.5) ** 2
+    return Constant(1.0) + alpha / cosh(beta * (r - gamma)) ** 2
+
+
+# With an initial mesh and a monitor function, we are able to construct a
+# :class:`~.MongeAmpereMover` instance and adapt the mesh.
+
+mover = MongeAmpereMover(mesh, ring_monitor, method="quasi_newton")
+mover.move()
+
+# The adapted mesh can be accessed via the `mesh` attribute of the mover. Plotting it,
+# we see that the adapted mesh has its resolution focused around a ring, as expected.
+
+fig, axes = plt.subplots()
+triplot(mover.mesh, axes=axes)
+axes.set_aspect(1)
+plt.savefig("monge_ampere1-adapted_mesh.jpg")
+
+# .. figure:: monge_ampere1-adapted_mesh.jpg
+#    :figwidth: 60%
+#    :align: center
+#
+# This tutorial can be dowloaded as a `Python script <monge_ampere1.py>`__.
+#
+#
+# .. rubric:: References
+#
+# .. bibliography:: demo_references.bib
+#    :filter: docname in docnames

movement/monge_ampere.py

@@ -392,9 +392,9 @@ def __init__(self, mesh, monitor_function, **kwargs):
         # Create functions to hold solution data
         self.V = self.P1*self.P1_ten
         self.phisigma = firedrake.Function(self.V)
-        self.phi, self.sigma = self.phisigma.split()
+        self.phi, self.sigma = self.phisigma.subfunctions
         self.phisigma_old = firedrake.Function(self.V)
-        self.phi_old, self.sigma_old = self.phisigma_old.split()
+        self.phi_old, self.sigma_old = self.phisigma_old.subfunctions
 
         # Initialise phi and sigma
         self.apply_initial_guess(**kwargs)

test/test_demos.py

@@ -0,0 +1,41 @@
+"""
+Checks that all demo scripts run.
+"""
+import glob
+import os
+from os.path import splitext
+import pytest
+import shutil
+import subprocess
+import sys
+
+
+# Set environment variable to specify that tests are being run so that we can
+# cut down the length of the demos
+os.environ["MOVEMENT_REGRESSION_TEST"] = "1"
+
+cwd = os.path.abspath(os.path.dirname(__file__))
+demo_dir = os.path.abspath(os.path.join(cwd, "..", "demos"))
+all_demos = glob.glob(os.path.join(demo_dir, "*.py"))
+
+
+@pytest.fixture(params=all_demos, ids=lambda x: splitext(x.split("demos/")[-1])[0])
+def demo_file(request):
+    return os.path.abspath(request.param)
+
+
+def test_demos(demo_file, tmpdir, monkeypatch):
+    assert os.path.isfile(demo_file), f"Demo file '{demo_file}' not found."
+
+    # Copy mesh files
+    source = os.path.dirname(demo_file)
+    for f in glob.glob(os.path.join(source, "*.msh")):
+        shutil.copy(f, str(tmpdir))
+
+    # Change working directory to temporary directory
+    monkeypatch.chdir(tmpdir)
+    subprocess.check_call([sys.executable, demo_file])
+
+    # Clean up plots
+    for ext in ("jpg", "pvd", "vtu"):
+        subprocess.check_call(["rm", "-f", f"*.{ext}"])

test/test_smoothing.py

@@ -63,7 +63,6 @@ def update_forcings(t):
         pwd = os.path.dirname(__file__)
         fname = os.path.join(pwd, "data", "forced_mesh_laplacian.npy")
         assert np.allclose(new_coords, np.load(fname))
-    return mover
 
 
 if __name__ == "__main__":

test/test_spring.py

@@ -52,7 +52,7 @@ def test_angles():
     assert np.allclose(angles, expected)
 
 
-def test_forced(method, time, plot=False, test=True):
+def test_forced(method, time, plot=False):
     """
     Test that a uniform mesh is moved as expected
     when only one of its boundaries is forced.
@@ -88,29 +88,22 @@ def update_forcings(t):
     # Plotting
     if plot:
         import matplotlib.pyplot as plt
+        from firedrake.pyplot import triplot
 
         fig, axes = plt.subplots()
-        firedrake.triplot(mover.mesh, axes=axes)
+        triplot(mover.mesh, axes=axes)
         axes.axis(False)
         plt.tight_layout()
         plt.savefig(f"plots/mesh_{method}_{it}.png")
 
     # Check as expected
-    if test:
-        pwd = os.path.dirname(__file__)
-        fname = os.path.join(pwd, "data", f"forced_mesh_lineal_{it}.npy")
-        expected = coords if np.isclose(time, 0.0) else np.load(fname)
-        assert np.allclose(new_coords, expected)
-    return mover
+    pwd = os.path.dirname(__file__)
+    fname = os.path.join(pwd, "data", f"forced_mesh_lineal_{it}.npy")
+    expected = coords if np.isclose(time, 0.0) else np.load(fname)
+    assert np.allclose(new_coords, expected)
 
-
-def test_plex_consistency(method, time):
-    """
-    Test that the Firedrake mesh coordinates and the
-    underlying DMPlex's coordinates are consistent
-    after mesh movement.
-    """
-    mover = test_forced(method, time, test=False)
+    # Test that the Firedrake mesh coordinates and the underlying DMPlex's
+    # coordinates are consistent after mesh movement.
     plex_coords = mover.plex.getCoordinatesLocal().array
     mesh_coords = mover.mesh.coordinates.dat.data_with_halos
     assert np.allclose(plex_coords.reshape(*mesh_coords.shape), mesh_coords)