From b7cde739ff9a8af247c3ac40f3f3cf716f8be22d Mon Sep 17 00:00:00 2001 From: Davor Dundovic <33790330+ddundo@users.noreply.github.com> Date: Thu, 13 Jun 2024 19:50:08 +0200 Subject: [PATCH] Minor demo and readme fixes (#98) I was going through movement and saw minor things that needed to be fixed: - broken link in README.md - typos and formatting issues in demos - cross-referencing in demos --- README.md | 4 ++-- demos/laplacian_smoothing.py | 14 +++++++++----- demos/lineal_spring.py | 15 +++++++++++++-- demos/monge_ampere1.py | 24 +++++++++++++----------- demos/monge_ampere_3d.py | 13 +++++++------ demos/monge_ampere_helmholtz.py | 16 ++++++++-------- 6 files changed, 52 insertions(+), 34 deletions(-) diff --git a/README.md b/README.md index 0ae765a..4079e14 100644 --- a/README.md +++ b/README.md @@ -10,9 +10,9 @@ Movement is built for use with the [Firedrake finite element library](https://fi ## Documentation Movement has two sources of documentation: -* The [website](mesh-adaptation.github.io), which includes long-form documentation, demos, and API documentation. +* The [website](https://mesh-adaptation.github.io/movement/index.html), which includes long-form documentation, demos, and API documentation. * The [wiki](https://github.com/mesh-adaptation/mesh-adaptation-docs/wiki), which includes recommendations on how to interact with the codebase and development practices to be followed. ## Installation -Forinstallation instructions, we refer to the [Wiki page](https://github.com/mesh-adaptation/mesh-adaptation-docs/wiki/Installation-Instructions). +For installation instructions, we refer to the [Wiki page](https://github.com/mesh-adaptation/mesh-adaptation-docs/wiki/Installation-Instructions). diff --git a/demos/laplacian_smoothing.py b/demos/laplacian_smoothing.py index 0e59cb2..339ff1f 100644 --- a/demos/laplacian_smoothing.py +++ b/demos/laplacian_smoothing.py @@ -32,7 +32,7 @@ # i.e., # # .. math:: -# \frac1{h^2}(-v_{i-1} + 2v_i - v_{i+1})) = 0, +# \frac1{h^2}(-v_{i-1} + 2v_i - v_{i+1}) = 0, # # the left-hand side of which you might recognise as a finite difference approximation # of the second derivative, i.e., the Laplace operator. @@ -61,13 +61,13 @@ # :align: center # # Suppose we wish to enforce a time-dependent velocity :math:`\mathbf{v}_f` -# on the top boundary :math and see how the mesh responds. Consider the velocity +# on the top boundary and see how the mesh responds. Consider the velocity # # .. math:: # \mathbf{v}_f(x,y,t) = \left[0, A\:\sin\left(\frac{2\pi t}T\right)\:\sin(\pi x)\right] # # acting only in the vertical direction, where :math:`A` is the amplitude and :math:`T` -# is the time period. The displacement following a sinusoidal pattern along that +# is the time period. The displacement follows a sinusoidal pattern along that # boundary. The boundary movement is also sinusoidal in time, such that it ramps up and # then reverses, with the analytical solution being that the final mesh coincides with # the initial mesh. @@ -103,8 +103,8 @@ def boundary_velocity(x, t): # # To enforce this boundary velocity, we need to create a :class:`~.LaplacianSmoother` # instance and define a function for updating the boundary conditions. Since we are -# going to enforce the velocity on the top boundary, we create a :class:`~.Function` to -# represent the boundary condition values and pass this to a :class:`~.DirichletBC` +# going to enforce the velocity on the top boundary, we create a :class:`~firedrake.function.Function` to +# represent the boundary condition values and pass this to a :class:`~firedrake.bcs.DirichletBC` # object. We then define a function which updates it as time progresses. :: mover = LaplacianSmoother(mesh, timestep) @@ -166,4 +166,8 @@ def update_boundary_velocity(t): print(f"l_infinity error: {linf_error:.3f} m") assert np.isclose(linf_error, 0.0) +# .. code-block:: console +# +# l_infinity error: 0.000 m +# # This tutorial can be downloaded as a `Python script `__. diff --git a/demos/lineal_spring.py b/demos/lineal_spring.py index 7c07947..0f7272a 100644 --- a/demos/lineal_spring.py +++ b/demos/lineal_spring.py @@ -122,8 +122,9 @@ def boundary_displacement(x, t): # :figwidth: 60% # :align: center # -# To apply this boundary displacement, we need to create a :class:`~.SpringMover` -# instance and define a function for updating the boundary conditions. :: +# To apply this boundary displacement, we need to create a +# :class:`~movement.spring.SpringMover` instance and define a function for updating the +# boundary conditions. :: mover = SpringMover(mesh, timestep, method="lineal") top = Function(mover.coord_space) @@ -182,6 +183,10 @@ def update_boundary_displacement(t): print(f"l_infinity error: {linf_error:.3f} m") assert linf_error < 0.002 +# .. code-block:: console +# +# l_infinity error: 0.001 m +# # Note that the mesh doesn't return to its original state quite as neatly with the lineal # spring method as it does with the Laplacian smoothing method. However, the result is # still very good (as can be seen from the plots above). @@ -192,6 +197,12 @@ def update_boundary_displacement(t): print(f"Stiffness matrix shape: {K.shape}") print(f"Number of mesh vertices: {mesh.num_vertices()}") +# .. code-block:: console +# +# Stiffness matrix shape: (242, 242) +# Number of mesh vertices: 121 +# + fig, axes = plt.subplots() axes.spy(K) plt.savefig("lineal_spring-sparsity.jpg") diff --git a/demos/monge_ampere1.py b/demos/monge_ampere1.py index b615e1f..1e40abd 100644 --- a/demos/monge_ampere1.py +++ b/demos/monge_ampere1.py @@ -19,14 +19,15 @@ # 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 +# the computational and physical domains, respectively. 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, +# such that it 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, @@ -34,20 +35,21 @@ # 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.) +# 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. +# We begin the example by importing from the namespaces of Firedrake and Movement. :: -# To start with a simple example, consider a uniform mesh of the unit square. Feel free -# to ignore the `"MOVEMENT_REGRESSION_TEST"`, as it is only used when this demo is run -# in the test suite (to reduce its runtime). :: import os from firedrake import * from movement import * +# To start with a simple example, consider a uniform mesh of the unit square. Feel free +# to ignore the `"MOVEMENT_REGRESSION_TEST"`, as it is only used when this demo is run +# in the test suite (to reduce its runtime). :: + test = os.environ.get("MOVEMENT_REGRESSION_TEST") n = 10 if test else 20 mesh = UnitSquareMesh(n, n) @@ -91,7 +93,7 @@ def ring_monitor(mesh): # With an initial mesh and a monitor function, we are able to construct a -# :class:`~.MongeAmpereMover` instance and adapt the mesh. By default, the Monge-Ampère +# :class:`~movement.monge_ampere.MongeAmpereMover` instance and adapt the mesh. By default, the Monge-Ampère # equation is solved to a relative tolerance of :math:`10^{-8}`. However, for the # purposes of continuous integration testing, a tolerance of :math:`10^{-3}` is used # instead to further reduce the runtime. :: @@ -134,7 +136,7 @@ def ring_monitor(mesh): # 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. +# that the Monge-Ampère method can also be applied in three dimensions. # # This tutorial can be dowloaded as a `Python script `__. # diff --git a/demos/monge_ampere_3d.py b/demos/monge_ampere_3d.py index 4ea4c5a..87e0433 100644 --- a/demos/monge_ampere_3d.py +++ b/demos/monge_ampere_3d.py @@ -15,7 +15,7 @@ def sinatan3(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 +# We will now try to use mesh movement to optimise 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 @@ -27,13 +27,14 @@ def sinatan3(mesh): # 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 +# of :math:`H`. We have normalised 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 +# means that 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 +# :class:`~movement.monge_ampere.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) @@ -92,8 +93,8 @@ def monitor(mesh): # :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. +# how to optimise the mesh for the discretisation of a PDE with the aim to +# minimise its discretisation error. # # This tutorial can be dowloaded as a `Python script `__. # diff --git a/demos/monge_ampere_helmholtz.py b/demos/monge_ampere_helmholtz.py index 3e1f063..83b3079 100644 --- a/demos/monge_ampere_helmholtz.py +++ b/demos/monge_ampere_helmholtz.py @@ -1,8 +1,8 @@ -# Using mesh movement to optimize the mesh for PDE solution +# Using mesh movement to optimise the mesh for PDE solution # ========================================================= # -# In this demo will demonstrate how we might use mesh movement to obtain a mesh -# that is optimized for solving a particular PDE. The general idea is that we +# In this demo we will demonstrate how we might use mesh movement to obtain a mesh +# that is optimised for solving a particular PDE. The general idea is that we # want to reduce numerical error by increasing resolution where the local error # is (expected to be) high and decrease it elsewhere. # @@ -27,12 +27,12 @@ # f(x, y) &= -\nabla^2 u(x, y) + u(x, y) # = \left[ -\frac{4(x-x_0)^2 + 4(y-y_0)^2}{w^4} # + \frac{4}{w^2} + 1 \right] -# \exp\left(-\frac{(x-x_0)^2 + (y-y_0)^2}{w^2}\right) +# \exp\left(-\frac{(x-x_0)^2 + (y-y_0)^2}{w^2}\right), # # where :math:`(x_0, y_0)` is the centre of the Gaussian with width :math:`w`. # Note that here we *first* choose the solution :math:`u` after which we can # compute what rhs :math:`f` should be, by substitution in the PDE, in order -# for :math:`u` to be the analytical solution. This so called Method of +# for :math:`u` to be the analytical solution. This so-called Method of # Manufactured Solutions approach is an easy way to construct PDE # configurations for which we know the analytical solution, e.g. for testing # purposes. @@ -101,7 +101,7 @@ 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 +# We will now try to use mesh movement to optimise the mesh to reduce # 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. @@ -167,7 +167,7 @@ def monitor(mesh): # 17 Min/Max 8.6775e-01 Residual 5.8585e-08 Variation (σ/μ) 9.9835e-01 # Converged in 17 iterations. # -# Plotting the resulting mesh +# Plotting the resulting mesh: from firedrake.pyplot import triplot @@ -193,7 +193,7 @@ def monitor(mesh): # # Of course, in many practical problems we do not actually have access # to the exact solution. We can then use the Hessian of the numerical -# solution in the monitor function. Calculating the Hessian we have to +# solution in the monitor function. When calculating the Hessian we have to # be a bit careful however: since our numerical FEM solution :math:`u_h` # is a piecewise linear function, its first gradient results in a # piecewise *constant* function, a vector-valued constant