From 398fe4836f052a0e377950260b201abe8e94ac3e Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 16:43:59 +0100 Subject: [PATCH 1/5] Plate bending demo using MITC Solve a Reissner-Mindlin plate with the Mixed Interpolation of Tensorial Components formulation, which avoids shear locking in the thin-plate limit by projecting the rotations into an edge-conforming Nedelec space. The reduction operator is written as a symbolic interpolate inside the variational form. --- demos/plate_bending_mitc/elements.svg | 325 ++++++++++++++++++ demos/plate_bending_mitc/elements.tex | 38 ++ .../plate_bending_mitc.py.rst | 172 +++++++++ docs/source/advanced_tut.rst | 1 + tests/firedrake/demos/test_demos_run.py | 1 + 5 files changed, 537 insertions(+) create mode 100644 demos/plate_bending_mitc/elements.svg create mode 100644 demos/plate_bending_mitc/elements.tex create mode 100644 demos/plate_bending_mitc/plate_bending_mitc.py.rst diff --git a/demos/plate_bending_mitc/elements.svg b/demos/plate_bending_mitc/elements.svg new file mode 100644 index 0000000000..618a7b0d67 --- /dev/null +++ b/demos/plate_bending_mitc/elements.svg @@ -0,0 +1,325 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + P1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + P1-iso-P2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Nedelec1 + + + + + + + + + + + + diff --git a/demos/plate_bending_mitc/elements.tex b/demos/plate_bending_mitc/elements.tex new file mode 100644 index 0000000000..915e47d0a5 --- /dev/null +++ b/demos/plate_bending_mitc/elements.tex @@ -0,0 +1,38 @@ +\documentclass[tikz,border=2pt]{standalone} +\usepackage[utf8]{inputenc} +\begin{document} +\begin{tikzpicture} + + % P1 Element - Filled + \begin{scope}[shift={(0,0)}] + \filldraw[fill=blue!20, draw=black, thick] (0,0) -- (2,0) -- (1,1.732) -- cycle; + \foreach \p in {(0,0), (2,0), (1,1.732)} \fill[black] \p circle (2pt); + \node at (1,-0.5) {P1}; + \end{scope} + + % P1-iso-P2 Element - Filled + \begin{scope}[shift={(4,0)}] + \filldraw[fill=green!20, draw=black, thick] (0,0) -- (2,0) -- (1,1.732) -- cycle; + \draw[draw=black, thick] (1,0) -- (1.5,0.866) -- (0.5,0.866) -- cycle; + \draw[draw=black, thick] (0,0) -- (0.5,0.866); + \draw[draw=black, thick] (2,0) -- (1.5,0.866); + \draw[draw=black, thick] (1,1.732) -- (0.5,0.866); + \draw[draw=black, thick] (1,1.732) -- (1.5,0.866); + \draw[draw=black, thick] (1,0) -- (0.5,0.866); + \draw[draw=black, thick] (1,0) -- (1.5,0.866); + \foreach \p in {(0,0), (2,0), (1,1.732), (1,0), (0.5,0.866), (1.5,0.866)} \fill[black] \p circle (2pt); + \node at (1,-0.5) {P1-iso-P2}; + \end{scope} + + % Nédélec (1st kind) Element - Filled + \begin{scope}[shift={(8,0)}] + \filldraw[fill=red!20, draw=black, thick] (0,0) -- (2,0) -- (1,1.732) -- cycle; + % Vectors centered along the edges, always black + \draw[->, ultra thick, black] (0.25, 0) -- (1.75, 0); + \draw[->, ultra thick, black] (1.75, 0.433) -- (1.25, 1.299); + \draw[->, ultra thick, black] (0.75, 1.299) -- (0.25, 0.433); + \node at (1,-0.5) {Nedelec 1}; + \end{scope} + +\end{tikzpicture} +\end{document} diff --git a/demos/plate_bending_mitc/plate_bending_mitc.py.rst b/demos/plate_bending_mitc/plate_bending_mitc.py.rst new file mode 100644 index 0000000000..ca7cd92b39 --- /dev/null +++ b/demos/plate_bending_mitc/plate_bending_mitc.py.rst @@ -0,0 +1,172 @@ +Plate Bending Using Mixed Interpolation of Tensorial Components (MITC) +====================================================================== + +This demo illustrates how to solve a Reissner-Mindlin plate problem using the +Mixed Interpolation of Tensorial Components (MITC) formulation. The main goal of +the MITC method is to prevent shear locking in the thin-plate limit (:math:`t \to 0`) +by projecting the rotation field into an edge-conforming Nédélec space. + +Background and Formulation +-------------------------- + +We model the plate using the Reissner-Mindlin equations. Given a domain +:math:`\Omega` representing the plate's mid-surface, we seek the transverse deflection +:math:`w` and the rotations :math:`\boldsymbol{\beta}` of the mid-surface normal. + +To avoid shear locking, the standard shear strain term :math:`(\nabla w - \boldsymbol{\beta})` +is modified. We project the rotation field :math:`\boldsymbol{\beta}` onto an +:math:`H(\text{curl})`-conforming space :math:`\boldsymbol{R}_h` (a Nédélec space) using a reduction +operator :math:`\Pi_h`. The shear strain term is then evaluated as: + +.. math:: + + \bar{\boldsymbol{\gamma}}_h = \nabla w - \Pi_h \boldsymbol{\beta} + +Element Spaces +-------------- + +The following diagram illustrates the element spaces used: + +.. image:: elements.svg + :align: center + +* **P1**: Standard linear Lagrange element for the deflection :math:`w`. +* **P1-iso-P2**: Constructed as macroelement, where the master triangle is divided into four smaller sub-triangles. This structure is utilized for the rotation field :math:`\boldsymbol{\beta}`. +* **Nédélec 1**: An :math:`H(\text{curl})`-conforming space used for the MITC projection. The blue arrows represent the edge-based degrees of freedom. + +Implementation +-------------- + +Thanks to native compilation support for symbolic interpolation nodes inside +variational forms, we can define the reduction operator directly using +Firedrake's symbolic ``interpolate`` function within the UFL expression. Under the +hood, this bypasses ``BaseFormAssembler`` specifically for the case where +we have a ``Form`` with an ``Interpolate`` node onto the same mesh. + +We begin by importing the Firedrake namespace. + +:: + + from firedrake import * + +Mesh and Geometry +----------------- + +First, we set up a simple mesh of a unit square representing our plate. + +:: + + n = 16 + mesh = UnitSquareMesh(n, n) + +Material and Physical Parameters +-------------------------------- + +We set up standard parameters for a thin plate. Here, we define the thickness +:math:`t`, Young's modulus :math:`E`, Poisson's ratio :math:`\nu`, and the shear correction factor +:math:`k_s`. + +:: + + t = Constant(0.01) + E = Constant(1e3) + nu = Constant(0.3) + k_s = Constant(5.0/6.0) + +The bending stiffness :math:`D` and shear stiffness :math:`G` are derived below. + +:: + + D = E * t**3 / (12 * (1 - nu**2)) + G = E / (2 * (1 + nu)) + G_shear = k_s * G * t + +Function Spaces +--------------- + +We construct a mixed function space for the deflection and the rotation. We use +linear Lagrange elements for the deflection and P1-iso-P2 elements for the rotation. +Crucially, we also define a Nédélec space :math:`R` of degree 1 to serve as our +target edge-conforming space for the MITC projection. Because the facets are +split by the P1-iso-P2 macro-triangulation, the Nédélec integral moments must +employ a composite quadrature scheme (``quad_scheme="iso"``). + +:: + + W = FunctionSpace(mesh, "Lagrange", 1) + B = VectorFunctionSpace(mesh, "Lagrange", 1, variant="iso") + R = FunctionSpace(mesh, "N1curl", 1, quad_scheme="iso") + + V = MixedFunctionSpace([W, B]) + +Variational Formulation +----------------------- + +Next, we declare the trial and test functions. We write the isotropic bending +stress tensor :math:`\boldsymbol{\sigma}(\boldsymbol{\beta})` and establish our bilinear forms. + +:: + + u = Function(V) + w, beta = TrialFunctions(V) + v, theta = TestFunctions(V) + + def sigma(phi): + return D * ((1 - nu) * sym(grad(phi)) + nu * div(phi) * Identity(2)) + + a_bending = inner(sigma(beta), sym(grad(theta))) * dx + +We implement the MITC projection using Firedrake's symbolic ``interpolate()`` +function. This acts as a true symbolic UFL operator embedded +directly within the form definition. TSFC handles the assembly +pipeline seamlessly by evaluating the node on the same mesh. + +:: + + Pi_beta = interpolate(beta, R) + Pi_theta = interpolate(theta, R) + + a_shear = G_shear * inner(grad(w) - Pi_beta, grad(v) - Pi_theta) * dx + + a = a_bending + a_shear + +We impose a uniform transverse downward load :math:`f` acting on the plate. + +:: + + f = Constant(1.0) + L = inner(f, v) * dx + +Boundary Conditions +------------------- + +We apply fully clamped boundary conditions on all boundaries, meaning both the +deflection :math:`w` and the rotation :math:`\boldsymbol{\beta}` vanish. + +:: + + bcs = [DirichletBC(V.sub(0), 0, "on_boundary"), + DirichletBC(V.sub(1), 0, "on_boundary")] + +Computation +----------- + +We solve the problem in the usual way. + +:: + + solve(a == L, u, bcs=bcs) + +We recover the split deflection and rotation solutions for analysis. + +:: + + w_sol, beta_sol = u.subfunctions + max_w = w_sol.dat.data.max() + print(f"Max deflection: {max_w:.6e}") + +Finally, we output the deflection to a PVD file for visualization in ParaView. + +:: + + VTKFile("mitc_plate.pvd").write(w_sol) diff --git a/docs/source/advanced_tut.rst b/docs/source/advanced_tut.rst index 5f3bc49bf6..0926f1178f 100644 --- a/docs/source/advanced_tut.rst +++ b/docs/source/advanced_tut.rst @@ -39,3 +39,4 @@ element systems. Nonlinear preconditioning using an auxiliary SNES for the Allen-Cahn equation. Reynolds-robust preconditioning of the stationary Navier-Stokes equations. Eigenvalue problem with guaranteed bounds and adaptive refinement. + A plate bending problem solved with Mixed Interpolation of Tensorial Components. diff --git a/tests/firedrake/demos/test_demos_run.py b/tests/firedrake/demos/test_demos_run.py index 236878a1e7..f8e29bd048 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -57,6 +57,7 @@ Demo(('submesh_reaction_diffusion', 'submesh_reaction_diffusion'), ["netgen", "vtk"]), Demo(('nonlinear_pc', 'nonlinear_pc_allen_cahn'), []), Demo(('reynolds_robust_navier_stokes_hdiv', 'reynolds_robust_navier_stokes_hdiv'), ["vtk"]), + Demo(('plate_bending_mitc', 'plate_bending_mitc'), ["vtk"]), ] PARALLEL_DEMOS = [ Demo(("full_waveform_inversion", "full_waveform_inversion"), ["adjoint"]), From bb3797bac7e583df68a17fc76b2262e05c1ef7a4 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 17:13:25 +0100 Subject: [PATCH 2/5] Modified Morley demo for a singular perturbation problem Solve the fourth order problem of Wang, Xu and Hu (2006) with the modified Morley element, which keeps Morley for the biharmonic term and replaces the Laplacian term by its linear conforming interpolant so that the method stays uniform as the perturbation parameter vanishes. The reduction operator is a symbolic interpolate of an argument sitting under a gradient inside the form, so it exercises a case the MITC demo does not: the form compiler has to differentiate the fused interpolant in place. The demo checks the degenerate limit, where the reduced solution is the linear conforming solution of the Poisson problem the equation approaches. Boundary conditions are imposed weakly because Firedrake has no strong boundary conditions on the Morley element. --- demos/modified_morley/modified_morley.py.rst | 177 +++++++++++++++++++ 1 file changed, 177 insertions(+) create mode 100644 demos/modified_morley/modified_morley.py.rst diff --git a/demos/modified_morley/modified_morley.py.rst b/demos/modified_morley/modified_morley.py.rst new file mode 100644 index 0000000000..51f21678fe --- /dev/null +++ b/demos/modified_morley/modified_morley.py.rst @@ -0,0 +1,177 @@ +Modified Morley Element for a Fourth-Order Singular Perturbation Problem +======================================================================= + +This demo solves a fourth-order elliptic singular perturbation problem with the +modified Morley element of Wang, Xu and Hu (*Modified Morley element method for +a fourth order elliptic singular perturbation problem*, J. Comput. Math. 24 +(2006), 113-120). Like the :doc:`MITC plate bending demo `, +it writes a reduction operator as a symbolic ``interpolate`` inside the +variational form; here the operator sits in the *lower* order term, and it is +what makes the method converge uniformly as the perturbation parameter vanishes. + +Background and Formulation +-------------------------- + +On a polygonal domain :math:`\Omega \subset \mathbb{R}^2` we seek :math:`u` with + +.. math:: + + \varepsilon^2 \Delta^2 u - \Delta u = f \quad \text{in } \Omega, \qquad + u = \frac{\partial u}{\partial \nu} = 0 \quad \text{on } \partial\Omega, + +for a small parameter :math:`0 < \varepsilon \le 1`. As :math:`\varepsilon \to 0` +the equation formally degenerates to the Poisson problem +:math:`-\Delta u^0 = f`, and a method for the fourth-order problem is only +useful here if it degenerates the same way. + +The Morley element is the cheapest triangular element for fourth-order +problems, but it is not a :math:`C^0` element: the plain Morley discretisation +of the equation above is divergent as :math:`\varepsilon \to 0`. The modified +method keeps the Morley element for the fourth-order term and replaces the +second-order term by its linear conforming interpolant. Writing +:math:`\Pi_h : V_h \to P_h` for interpolation into the linear Lagrange space, +the discrete problem is + +.. math:: + + \varepsilon^2 a_h(u_h, v_h) + b_h(\Pi_h u_h, \Pi_h v_h) + = (f, \Pi_h v_h) \qquad \forall v_h \in V_h, + +with the broken forms + +.. math:: + + a_h(v, w) = \sum_{T} \int_T \nabla^2 v : \nabla^2 w \, \mathrm{d}x, \qquad + b_h(v, w) = \sum_{T} \int_T \nabla v \cdot \nabla w \, \mathrm{d}x. + +At :math:`\varepsilon = 0` this collapses to :math:`b_h(\Pi_h u_h, \Pi_h v_h) = +(f, \Pi_h v_h)`, which is the linear conforming discretisation of the Poisson +problem -- so :math:`\Pi_h u_h` is exactly the :math:`P_1` solution of the +degenerate equation. We reproduce that limit numerically at the end. + +We begin by importing the Firedrake namespace. + +:: + + from firedrake import * + +Mesh and Function Spaces +------------------------ + +We take a uniform triangulation of the unit square. :math:`V` carries the +Morley element and :math:`P` the linear conforming space that the reduction +operator maps into. + +:: + + n = 16 + mesh = UnitSquareMesh(n, n) + + V = FunctionSpace(mesh, "Morley", 2) + P = FunctionSpace(mesh, "CG", 1) + +Variational Formulation +----------------------- + +We declare the trial and test functions on the Morley space, along with the +perturbation parameter. + +:: + + u = TrialFunction(V) + v = TestFunction(V) + + epsilon = Constant(1E-2) + +The reduction operator is written directly as a symbolic ``interpolate`` of an +argument. Because the interpolation is onto the same mesh, the form compiler +fuses it into the surrounding integral rather than assembling a separate +operator, so ``grad(Pi_u)`` differentiates the interpolant in place. + +:: + + Pi_u = interpolate(u, P) + Pi_v = interpolate(v, P) + + a_bending = inner(grad(grad(u)), grad(grad(v))) * dx + a_membrane = inner(grad(Pi_u), grad(Pi_v)) * dx + +Boundary Conditions +------------------- + +The clamped conditions ask for both :math:`u` and its normal derivative to +vanish. Firedrake does not implement strong boundary conditions on the Morley +element, whose degrees of freedom are not all point values, so we impose both +weakly with a penalty. The deflection is penalised through the conforming +interpolant, which is the trace the second-order term sees, and the normal +derivative is penalised on the Morley function itself. + +:: + + normal = FacetNormal(mesh) + h = CellDiameter(mesh) + alpha = Constant(20.0) + + a_penalty = (alpha / h * inner(Pi_u, Pi_v) * ds + + epsilon**2 * alpha / h + * inner(dot(grad(u), normal), dot(grad(v), normal)) * ds) + + a = epsilon**2 * a_bending + a_membrane + a_penalty + +The load is applied through the same reduction operator as the second-order +term, matching the right-hand side of the discrete problem. + +:: + + f = Constant(1.0) + L = inner(f, Pi_v) * dx + +Computation +----------- + +We solve the problem in the usual way. + +:: + + uh = Function(V) + solve(a == L, uh) + +The Degenerate Limit +-------------------- + +To see that the method degenerates correctly, we solve the Poisson problem that +the equation approaches, discretised with the same linear conforming space and +the same weak boundary condition. + +:: + + p = TrialFunction(P) + q = TestFunction(P) + + a_poisson = inner(grad(p), grad(q)) * dx + alpha / h * inner(p, q) * ds + L_poisson = inner(f, q) * dx + + u_poisson = Function(P) + solve(a_poisson == L_poisson, u_poisson) + +The reduced solution :math:`\Pi_h u_h` should agree with it to +:math:`O(\varepsilon^2)`. + +:: + + Pi_uh = assemble(interpolate(uh, P)) + difference = errornorm(u_poisson, Pi_uh) + print(f"epsilon = {float(epsilon):.1e}, " + f"relative difference = {difference / norm(u_poisson):.4f}") + + assert difference < 0.05 * norm(u_poisson) + +Repeating the solve over a range of :math:`\varepsilon` shows the difference +falling quadratically, so the discretisation is uniform in the perturbation +parameter rather than degenerating with it. + +Finally, we output the reduced deflection for visualisation in ParaView. + +:: + + VTKFile("modified_morley.pvd").write(Pi_uh) From 60ee5de5d2bb0f2b07d73e46048558c853c4508a Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 17:14:39 +0100 Subject: [PATCH 3/5] Draw the Morley and P1 elements in the modified Morley demo The reduction operator is easier to read as a picture: Morley carries vertex values and averaged edge normal derivatives, P1 carries the vertex values alone, and the operator keeps the ones they share. --- demos/modified_morley/elements.svg | 77 ++++++++++++++++++++ demos/modified_morley/elements.tex | 32 ++++++++ demos/modified_morley/modified_morley.py.rst | 15 ++++ docs/source/advanced_tut.rst | 1 + tests/firedrake/demos/test_demos_run.py | 1 + 5 files changed, 126 insertions(+) create mode 100644 demos/modified_morley/elements.svg create mode 100644 demos/modified_morley/elements.tex diff --git a/demos/modified_morley/elements.svg b/demos/modified_morley/elements.svg new file mode 100644 index 0000000000..d54301dc64 --- /dev/null +++ b/demos/modified_morley/elements.svg @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/demos/modified_morley/elements.tex b/demos/modified_morley/elements.tex new file mode 100644 index 0000000000..5f99eba909 --- /dev/null +++ b/demos/modified_morley/elements.tex @@ -0,0 +1,32 @@ +\documentclass[tikz,border=2pt]{standalone} +\usepackage[utf8]{inputenc} +\begin{document} +\begin{tikzpicture} + + % Morley Element - vertex values and edge normal derivatives + \begin{scope}[shift={(0,0)}] + \filldraw[fill=blue!20, draw=black, thick] (0,0) -- (2,0) -- (1,1.732) -- cycle; + % Outward normal derivative averages at the edge midpoints + \draw[->, ultra thick, black] (1,0) -- (1,-0.55); + \draw[->, ultra thick, black] (1.5,0.866) -- (1.976,1.141); + \draw[->, ultra thick, black] (0.5,0.866) -- (0.024,1.141); + % Point values at the vertices + \foreach \p in {(0,0), (2,0), (1,1.732)} \fill[black] \p circle (2pt); + \node at (1,-1.1) {Morley}; + \end{scope} + + % Reduction operator + \begin{scope}[shift={(0,0)}] + \draw[->, very thick, black] (2.6,0.866) -- (3.9,0.866); + \node at (3.25,1.25) {$\Pi_h$}; + \end{scope} + + % P1 Element - vertex values + \begin{scope}[shift={(4.5,0)}] + \filldraw[fill=green!20, draw=black, thick] (0,0) -- (2,0) -- (1,1.732) -- cycle; + \foreach \p in {(0,0), (2,0), (1,1.732)} \fill[black] \p circle (2pt); + \node at (1,-1.1) {P1}; + \end{scope} + +\end{tikzpicture} +\end{document} diff --git a/demos/modified_morley/modified_morley.py.rst b/demos/modified_morley/modified_morley.py.rst index 51f21678fe..bc2dec0613 100644 --- a/demos/modified_morley/modified_morley.py.rst +++ b/demos/modified_morley/modified_morley.py.rst @@ -49,6 +49,21 @@ At :math:`\varepsilon = 0` this collapses to :math:`b_h(\Pi_h u_h, \Pi_h v_h) = problem -- so :math:`\Pi_h u_h` is exactly the :math:`P_1` solution of the degenerate equation. We reproduce that limit numerically at the end. +Element Spaces +-------------- + +The following diagram shows the two elements and the operator between them: + +.. image:: elements.svg + :align: center + +* **Morley**: the quadratic nonconforming element for fourth-order problems. + Its degrees of freedom are the values at the vertices and the averaged + outward normal derivatives on the edges, drawn as arrows. +* **P1**: the linear Lagrange element, whose degrees of freedom are the vertex + values alone. The reduction operator :math:`\Pi_h` keeps those vertex values + and discards the normal derivatives. + We begin by importing the Firedrake namespace. :: diff --git a/docs/source/advanced_tut.rst b/docs/source/advanced_tut.rst index 0926f1178f..604f533180 100644 --- a/docs/source/advanced_tut.rst +++ b/docs/source/advanced_tut.rst @@ -40,3 +40,4 @@ element systems. Reynolds-robust preconditioning of the stationary Navier-Stokes equations. Eigenvalue problem with guaranteed bounds and adaptive refinement. A plate bending problem solved with Mixed Interpolation of Tensorial Components. + A fourth order singular perturbation problem solved with the modified Morley element. diff --git a/tests/firedrake/demos/test_demos_run.py b/tests/firedrake/demos/test_demos_run.py index f8e29bd048..2f661b6ae1 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -58,6 +58,7 @@ Demo(('nonlinear_pc', 'nonlinear_pc_allen_cahn'), []), Demo(('reynolds_robust_navier_stokes_hdiv', 'reynolds_robust_navier_stokes_hdiv'), ["vtk"]), Demo(('plate_bending_mitc', 'plate_bending_mitc'), ["vtk"]), + Demo(('modified_morley', 'modified_morley'), ["vtk"]), ] PARALLEL_DEMOS = [ Demo(("full_waveform_inversion", "full_waveform_inversion"), ["adjoint"]), From e7e3504e75fc02bafe65b4256c341c8f7d897af9 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 17:16:55 +0100 Subject: [PATCH 4/5] Give each demo's element diagram a distinct file name The docs Makefile copies every demo's assets into one flat directory, so two demos cannot both ship an elements.svg: the copy refuses to overwrite and the build stops before Sphinx runs. --- demos/modified_morley/modified_morley.py.rst | 4 ++-- demos/modified_morley/{elements.svg => morley_elements.svg} | 0 demos/modified_morley/{elements.tex => morley_elements.tex} | 0 demos/plate_bending_mitc/{elements.svg => mitc_elements.svg} | 0 demos/plate_bending_mitc/{elements.tex => mitc_elements.tex} | 0 demos/plate_bending_mitc/plate_bending_mitc.py.rst | 2 +- 6 files changed, 3 insertions(+), 3 deletions(-) rename demos/modified_morley/{elements.svg => morley_elements.svg} (100%) rename demos/modified_morley/{elements.tex => morley_elements.tex} (100%) rename demos/plate_bending_mitc/{elements.svg => mitc_elements.svg} (100%) rename demos/plate_bending_mitc/{elements.tex => mitc_elements.tex} (100%) diff --git a/demos/modified_morley/modified_morley.py.rst b/demos/modified_morley/modified_morley.py.rst index bc2dec0613..ca1e768deb 100644 --- a/demos/modified_morley/modified_morley.py.rst +++ b/demos/modified_morley/modified_morley.py.rst @@ -1,5 +1,5 @@ Modified Morley Element for a Fourth-Order Singular Perturbation Problem -======================================================================= +======================================================================== This demo solves a fourth-order elliptic singular perturbation problem with the modified Morley element of Wang, Xu and Hu (*Modified Morley element method for @@ -54,7 +54,7 @@ Element Spaces The following diagram shows the two elements and the operator between them: -.. image:: elements.svg +.. image:: morley_elements.svg :align: center * **Morley**: the quadratic nonconforming element for fourth-order problems. diff --git a/demos/modified_morley/elements.svg b/demos/modified_morley/morley_elements.svg similarity index 100% rename from demos/modified_morley/elements.svg rename to demos/modified_morley/morley_elements.svg diff --git a/demos/modified_morley/elements.tex b/demos/modified_morley/morley_elements.tex similarity index 100% rename from demos/modified_morley/elements.tex rename to demos/modified_morley/morley_elements.tex diff --git a/demos/plate_bending_mitc/elements.svg b/demos/plate_bending_mitc/mitc_elements.svg similarity index 100% rename from demos/plate_bending_mitc/elements.svg rename to demos/plate_bending_mitc/mitc_elements.svg diff --git a/demos/plate_bending_mitc/elements.tex b/demos/plate_bending_mitc/mitc_elements.tex similarity index 100% rename from demos/plate_bending_mitc/elements.tex rename to demos/plate_bending_mitc/mitc_elements.tex diff --git a/demos/plate_bending_mitc/plate_bending_mitc.py.rst b/demos/plate_bending_mitc/plate_bending_mitc.py.rst index ca7cd92b39..ba7b21ab82 100644 --- a/demos/plate_bending_mitc/plate_bending_mitc.py.rst +++ b/demos/plate_bending_mitc/plate_bending_mitc.py.rst @@ -27,7 +27,7 @@ Element Spaces The following diagram illustrates the element spaces used: -.. image:: elements.svg +.. image:: mitc_elements.svg :align: center * **P1**: Standard linear Lagrange element for the deflection :math:`w`. From a664c8d78f87f61a32d99a2f510dc1c69a9982d2 Mon Sep 17 00:00:00 2001 From: Pablo Brubeck Date: Thu, 27 Aug 2026 17:29:37 +0100 Subject: [PATCH 5/5] Cite the sources and link the runnable script Both demos rendered without the link to the script pylit generates from them, which is how a reader gets a version they can run. The modified Morley demo named its source in prose rather than citing it, so the references it rests on were absent from the bibliography. --- demos/demo_references.bib | 22 +++++++++++++++++++ demos/modified_morley/modified_morley.py.rst | 16 +++++++++----- .../plate_bending_mitc.py.rst | 2 ++ 3 files changed, 35 insertions(+), 5 deletions(-) diff --git a/demos/demo_references.bib b/demos/demo_references.bib index 4e1fe38dc7..7eabe788b0 100644 --- a/demos/demo_references.bib +++ b/demos/demo_references.bib @@ -601,3 +601,25 @@ @article{Hu:2014 doi={10.1007/s10915-014-9821-5} } + +@article{Wang:2006, + title={Modified {M}orley element method for a fourth order elliptic singular perturbation problem}, + author={Wang, Ming and Xu, Jin-chao and Hu, Yu-cheng}, + journal={Journal of Computational Mathematics}, + volume={24}, + number={2}, + pages={113--120}, + year={2006}, + url={https://www.jstor.org/stable/43694071} +} + +@article{Morley:1968, + title={The triangular equilibrium element in the solution of plate bending problems}, + author={Morley, L. S. D.}, + journal={Aeronautical Quarterly}, + volume={19}, + number={2}, + pages={149--169}, + year={1968}, + doi={10.1017/S0001925900004546} +} diff --git a/demos/modified_morley/modified_morley.py.rst b/demos/modified_morley/modified_morley.py.rst index ca1e768deb..643bb98fc5 100644 --- a/demos/modified_morley/modified_morley.py.rst +++ b/demos/modified_morley/modified_morley.py.rst @@ -2,9 +2,8 @@ Modified Morley Element for a Fourth-Order Singular Perturbation Problem ======================================================================== This demo solves a fourth-order elliptic singular perturbation problem with the -modified Morley element of Wang, Xu and Hu (*Modified Morley element method for -a fourth order elliptic singular perturbation problem*, J. Comput. Math. 24 -(2006), 113-120). Like the :doc:`MITC plate bending demo `, +modified Morley element of Wang, Xu and Hu :cite:`Wang:2006`. Like the +:doc:`MITC plate bending demo `, it writes a reduction operator as a symbolic ``interpolate`` inside the variational form; here the operator sits in the *lower* order term, and it is what makes the method converge uniformly as the perturbation parameter vanishes. @@ -24,8 +23,8 @@ the equation formally degenerates to the Poisson problem :math:`-\Delta u^0 = f`, and a method for the fourth-order problem is only useful here if it degenerates the same way. -The Morley element is the cheapest triangular element for fourth-order -problems, but it is not a :math:`C^0` element: the plain Morley discretisation +The Morley element :cite:`Morley:1968` is the cheapest triangular element for +fourth-order problems, but it is not a :math:`C^0` element: the plain Morley discretisation of the equation above is divergent as :math:`\varepsilon \to 0`. The modified method keeps the Morley element for the fourth-order term and replaces the second-order term by its linear conforming interpolant. Writing @@ -190,3 +189,10 @@ Finally, we output the reduced deflection for visualisation in ParaView. :: VTKFile("modified_morley.pvd").write(Pi_uh) + +A python script version of this demo can be found :demo:`here `. + +.. rubric:: References + +.. bibliography:: demo_references.bib + :filter: docname in docnames diff --git a/demos/plate_bending_mitc/plate_bending_mitc.py.rst b/demos/plate_bending_mitc/plate_bending_mitc.py.rst index ba7b21ab82..2539898675 100644 --- a/demos/plate_bending_mitc/plate_bending_mitc.py.rst +++ b/demos/plate_bending_mitc/plate_bending_mitc.py.rst @@ -170,3 +170,5 @@ Finally, we output the deflection to a PVD file for visualization in ParaView. :: VTKFile("mitc_plate.pvd").write(w_sol) + +A python script version of this demo can be found :demo:`here `.