diff --git a/.flake8 b/.flake8 new file mode 100644 index 0000000..9e59c6d --- /dev/null +++ b/.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 diff --git a/.gitignore b/.gitignore index d5750df..937bfea 100644 --- a/.gitignore +++ b/.gitignore @@ -1,7 +1,6 @@ -plots/ __pycache__/ movement.egg-info/ -*.ipynb_checkpoints/ +*.jpg *.pvd *.vtu diff --git a/Makefile b/Makefile index d42541c..31287a8 100644 --- a/Makefile +++ b/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 diff --git a/demos/demo_references.bib b/demos/demo_references.bib new file mode 100644 index 0000000..ddc7b41 --- /dev/null +++ b/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} +} diff --git a/demos/monge_ampere1.py b/demos/monge_ampere1.py new file mode 100644 index 0000000..d3fd2e3 --- /dev/null +++ b/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 `__. +# +# +# .. rubric:: References +# +# .. bibliography:: demo_references.bib +# :filter: docname in docnames diff --git a/movement/monge_ampere.py b/movement/monge_ampere.py index 7d6d28b..0bb6007 100644 --- a/movement/monge_ampere.py +++ b/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) diff --git a/test/test_demos.py b/test/test_demos.py new file mode 100644 index 0000000..c49378c --- /dev/null +++ b/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}"]) diff --git a/test/test_smoothing.py b/test/test_smoothing.py index 3ccebb0..aca0452 100644 --- a/test/test_smoothing.py +++ b/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__": diff --git a/test/test_spring.py b/test/test_spring.py index 81ede34..fb6d13a 100644 --- a/test/test_spring.py +++ b/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)