#58: Rework and reorder demos

demos/laplacian_smoothing.py

@@ -1,52 +1,58 @@
 # Mesh movement using Laplacian smoothing
 # =======================================
 #
-# In this demo, we again consider the simple boundary forcing problem from the
-# `lineal spring demo <lineal_spring.py.html>`__. However, instead of using the lineal
-# spring method, we use Laplacian smoothing.
-#
-# Recall that the lineal spring approach interprets the mesh as a fictitious structure
-# of beams, assigning the beams lengths and stiffness values. Given some forcing of the
-# beams, we determine how the rest of the structure responds by solving a discrete
-# Poisson equation. The difference with the Laplacian smoothing approach is that we do
-# not make use of the beam structure and simply formulate the Poisson equation in terms
-# of updates to the mesh coordinates, via a *mesh velocity*. That is, we solve
+# In this demo, we demonstrate the *Laplacian smoothing* approach. This method relies on
+# there being a user-specified forcing function, :math:`f`. With this, we can define a
+# Poisson equation of the form
 #
 # .. math::
 #    -\Delta\mathbf{v} = f
 #
-# and then compute
+# where :math:`\mathbf{v}` is the so-called *mesh velocity* that we solve for. That is,
+# with the mesh velocity, we update the mesh coordinates according to
 #
 # .. math::
-#    \mathbf{u} := \mathbf{u} + \mathbf{v}
+#    \mathbf{u} := \mathbf{u} + \mathbf{v} * \Delta t,
 #
-# where :math:`\mathbf{u}` is the mesh coordinate field, :math:`\mathbf{v}` is the mesh
-# velocity that we solve for under the forcing :math:`f`.
+# where :math:`\mathbf{u}` is the mesh coordinate field and :math:`\Delta t` is the
+# timestep.
 #
-# As ever, we begin by importing from the namespaces of Firedrake and Movement. We also
-# import various plotting utilites. ::
+# We begin by importing from the namespaces of Firedrake and Movement. ::
 
 from firedrake import *
 from movement import *
-from firedrake.pyplot import triplot
-import matplotlib.pyplot as plt
-import matplotlib.patches as mpatches
 
-# Recall the initial uniform mesh from the lineal spring demo, which is used again here.
-# ::
+# Let's start with a uniform mesh of the unit square. It has four boundary segments,
+# which are tagged with the integers 1, 2, 3, and 4. Note that segment 4 corresponds to
+# the top boundary. ::
+
+import matplotlib.pyplot as plt
+from firedrake.pyplot import triplot
 
-mesh = UnitSquareMesh(10, 10)
+n = 10
+mesh = UnitSquareMesh(n, n)
 coord_data_init = mesh.coordinates.dat.data.copy()
+fig, axes = plt.subplots()
+triplot(mesh, axes=axes)
+axes.set_aspect(1)
+axes.legend()
+plt.savefig("laplacian_smoothing-initial_mesh.jpg")
 
-# .. figure:: lineal_spring-initial_mesh.jpg
+# .. figure:: laplacian_smoothing-initial_mesh.jpg
 #    :figwidth: 50%
 #    :align: center
 #
-# In the lineal spring demo, the forcing corresponds to vertical
-# movement of the top boundary, with the displacement following a sinusoidal pattern
-# along that boundary. The forcing 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.
+# Suppose we wish to apply a time-dependent forcing to the top boundary and see how the
+# mesh responds. Consider the forcing
+#
+# .. math::
+#     \mathbf{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
+# boundary. The forcing 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.
 
 import numpy as np
 
@@ -60,17 +66,27 @@ def forcing(x, t):
     return forcing_amplitude * np.sin(2 * pi * t / forcing_period) * np.sin(pi * x)
 
 
+X = np.linspace(0, 1, n + 1)
 times = np.arange(0, forcing_period + 0.5 * timestep, timestep)
 
-# .. figure:: lineal_spring-forcings.jpg
+fig, axes = plt.subplots()
+for t in times:
+    axes.plot(X, forcing(X, t), label=f"t={t:.1f}")
+axes.set_xlim([0, 1])
+axes.legend()
+box = axes.get_position()
+axes.set_position([box.x0, box.y0, box.width * 0.8, box.height])
+axes.legend(loc="center left", bbox_to_anchor=(1, 0.5))
+plt.savefig("laplacian_smoothing-forcings.jpg")
+
+# .. figure:: laplacian_smoothing-forcings.jpg
 #    :figwidth: 60%
 #    :align: center
 #
-# The setup for the :class:`~.LaplacianSmoother` class is very similar to that for the
-# :class:`~.SpringMover class. The only difference is that we need to specify the
-# timestep value in addition to the mesh. Since we are going to apply a forcing to the
-# top boundary, we need to extract the indices for the associated boundary nodes. We
-# can then define the :func:`update_forcings` as in the other demo. ::
+# To apply this forcing, we need to create a :class:`~.LaplacianSmoother` instance and
+# define a function for updating the forcing applied to the boundary nodes. Since we are
+# going to apply a forcing to the top boundary, we need to extract the indices for the
+# associated boundary nodes. This is done using a :class:`~.DirichletBC` object. ::
 
 mover = LaplacianSmoother(mesh, timestep)
 boundary_nodes = DirichletBC(mesh.coordinates.function_space(), 0, 4).nodes
@@ -84,14 +100,18 @@ def update_forcings(t):
         forcing_data[i][1] = forcing(x, t)
 
 
-# We are now able to apply the mesh movement method. This works just as before. ::
+# We are now able to apply the mesh movement method. The forcings effectively enforce a
+# Dirichlet condition on the top boundary. On other boundaries, we enforce that there is
+# no movement using the `fixed_boundaries` keyword argument. ::
+
+import matplotlib.patches as mpatches
 
 time = 0.0
 fig, axes = plt.subplots(ncols=4, nrows=3, figsize=(12, 10))
 for i, t in enumerate(times):
     idx = 0 if i == 0 else (i + 1)
 
-    # Move the mesh and calculate the mesh speed
+    # Move the mesh and calculate the displacement
     mover.move(t, update_forcings=update_forcings, fixed_boundaries=[1, 2, 3])
     displacement = np.linalg.norm(mover.displacement)
     print(f"time = {t:.1f} s, displacement = {displacement:.2f} m")
@@ -108,9 +128,9 @@ def update_forcings(t):
 #    :figwidth: 100%
 #    :align: center
 #
-# Again, the mesh is deformed according to the vertical forcing, with the left, right,
-# and bottom boundaries remaining fixed, returning it to be very close to its original
-# state after one period. Let's check this in the :math:`\ell_\infty` norm. ::
+# The mesh is deformed according to the vertical forcing, with the left, right, and
+# bottom boundaries remaining fixed, returning it to be very close to its original state
+# after one period. Let's check this in the :math:`\ell_\infty` norm. ::
 
 coord_data = mover.mesh.coordinates.dat.data
 linf_error = np.max(np.abs(coord_data - coord_data_init))

demos/lineal_spring.py

@@ -2,7 +2,9 @@
 # ==============================================
 #
 # In this demo, we demonstrate a basic example using the *lineal spring* method, as
-# described in :cite:`Farhat:1998`.
+# described in :cite:`Farhat:1998`. For simplicity of presentation, we consider the same
+# example considered in the `Laplacian smoothing <laplacian_smoothing.py.html>`__ demo,
+# where mesh movement is driven by a forcing on the top boundary of a square mesh.
 #
 # The idea of the lineal spring method is to re-interpret the edges of a mesh as a
 # structure of stiff beams. Each beam has a stiffness associated with it, which is
@@ -22,7 +24,9 @@
 # :math:`\underline{\mathbf{K}_{ij}}\in\mathbb{R}^{2\times2}` and
 # :math:`\underline{\mathbf{K}}\in\mathbb{R}^{2N\times2N}`.
 #
-# Suppose we apply a forcing, which acts on the vertices according to a forcing matrix,
+# As with the Laplacian smoothing method, the lineal spring approach relies on there
+# being a user-specified forcing function. This acts on the vertices according to a
+# forcing matrix,
 #
 # .. math::
 #     \underline{\mathbf{f}} = \begin{bmatrix}
@@ -47,9 +51,9 @@
 from firedrake import *
 from movement import *
 
-# Let's start with a uniform mesh of the unit square. It has four boundary segments,
-# which are tagged with the integers 1, 2, 3, and 4. Note that segment 4 corresponds to
-# the top boundary. ::
+# Recall the initial uniform mesh of the unit square used in the Laplacian smoothing
+# demo, which has four boundary segments tagged with the integers 1, 2, 3, and 4. Note
+# that segment 4 corresponds to the top boundary. ::
 
 import matplotlib.pyplot as plt
 from firedrake.pyplot import triplot
@@ -67,14 +71,22 @@
 #    :figwidth: 50%
 #    :align: center
 #
-# Suppose we wish to apply a time-dependent forcing to the top boundary and see how the
-# mesh structure responds. Consider the forcing
+# We consider the same time-dependent forcing to the top boundary and see how the mesh
+# structure responds. Recall the formula
 #
 # .. math::
-#     \mathbf{f}(x,y,t)=\left[0, A\:\sin\left(\frac{2\pi t}T\right)\:\sin(\pi x)\right]
+#     \mathbf{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. ::
+# where :math:`A` is the amplitude and :math:`T` is the time period. However, we cannot
+# simply apply this formula in the same way as when solving the Poisson equation. We
+# first need to apply a transformation to account for the fact that we are considering
+# a discrete Poisson formulation. This amounts to applying the scale factor
+#
+# .. math::
+#    K = \Delta x^2 / 4,
+#
+# where :math:`\Delta x` is the local mesh size. To account for this, we need to provide
+# an additional index argument to :func:`forcing`. ::
 
 import numpy as np
 
@@ -108,8 +120,7 @@ def forcing(index, x, t):
 #    :align: center
 #
 # To apply this forcing, we need to create a :class:`~.SpringMover` instance and define
-# a function for updating the forcing applied to the boundary nodes. The way we get the
-# right indices for the top boundary is using a :class:`~.DirichletBC` object. ::
+# a function for updating the forcing applied to the boundary nodes. ::
 
 mover = SpringMover(mesh, method="lineal")
 
@@ -149,16 +160,19 @@ def update_forcings(t):
 #    :figwidth: 100%
 #    :align: center
 #
-# The mesh is deformed according to the vertical forcing, with the left, right, and
-# bottom boundaries remaining fixed, returning to be very close to its original state
+# Again, the mesh is deformed according to the vertical forcing, with the left, right,
+# and bottom boundaries remaining fixed, returning to be very close to its original state
 # after one period. Let's check this in the :math:`\ell_\infty` norm. ::
 
 coord_data = mover.mesh.coordinates.dat.data
 linf_error = np.max(np.abs(coord_data - coord_data_init))
 print(f"l_infinity error: {linf_error:.3f} m")
 assert linf_error < 0.01
 
-# Note that we can view the sparsity pattern of the stiffness matrix as follows. ::
+# 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.
+#
+# We can view the sparsity pattern of the stiffness matrix as follows. ::
 
 K = mover.stiffness_matrix
 print(f"Stiffness matrix shape: {K.shape}")