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 new file mode 100644 index 0000000000..643bb98fc5 --- /dev/null +++ b/demos/modified_morley/modified_morley.py.rst @@ -0,0 +1,198 @@ +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 :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. + +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 :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 +: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. + +Element Spaces +-------------- + +The following diagram shows the two elements and the operator between them: + +.. image:: morley_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. + +:: + + 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) + +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/modified_morley/morley_elements.svg b/demos/modified_morley/morley_elements.svg new file mode 100644 index 0000000000..d54301dc64 --- /dev/null +++ b/demos/modified_morley/morley_elements.svg @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/demos/modified_morley/morley_elements.tex b/demos/modified_morley/morley_elements.tex new file mode 100644 index 0000000000..5f99eba909 --- /dev/null +++ b/demos/modified_morley/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/plate_bending_mitc/mitc_elements.svg b/demos/plate_bending_mitc/mitc_elements.svg new file mode 100644 index 0000000000..618a7b0d67 --- /dev/null +++ b/demos/plate_bending_mitc/mitc_elements.svg @@ -0,0 +1,325 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + P1 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + P1-iso-P2 + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + Nedelec1 + + + + + + + + + + + + diff --git a/demos/plate_bending_mitc/mitc_elements.tex b/demos/plate_bending_mitc/mitc_elements.tex new file mode 100644 index 0000000000..915e47d0a5 --- /dev/null +++ b/demos/plate_bending_mitc/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..2539898675 --- /dev/null +++ b/demos/plate_bending_mitc/plate_bending_mitc.py.rst @@ -0,0 +1,174 @@ +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:: mitc_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) + +A python script version of this demo can be found :demo:`here `. diff --git a/docs/source/advanced_tut.rst b/docs/source/advanced_tut.rst index 5f3bc49bf6..604f533180 100644 --- a/docs/source/advanced_tut.rst +++ b/docs/source/advanced_tut.rst @@ -39,3 +39,5 @@ 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. + 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 236878a1e7..2f661b6ae1 100644 --- a/tests/firedrake/demos/test_demos_run.py +++ b/tests/firedrake/demos/test_demos_run.py @@ -57,6 +57,8 @@ 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"]), + Demo(('modified_morley', 'modified_morley'), ["vtk"]), ] PARALLEL_DEMOS = [ Demo(("full_waveform_inversion", "full_waveform_inversion"), ["adjoint"]),