Biharmonic equation

Authors: Julius Herb, Ottar Hellan and Jørgen S. Dokken

Download sources

• Python script

• Jupyter notebook

This demo illustrates how to:

• Solve a linear partial differential equation

• Use a discontinuous Galerkin method

• Solve a fourth-order differential equation

• Use the method of manufactured solutions

Equation and problem definition

Strong formulation

The biharmonic equation is a fourth-order elliptic equation. On the domain \(\Omega \subset \mathbb{R}^{d}\), \(1 \le d \le 3\), it reads

\[ \nabla^{4} u = f \quad {\rm in} \ \Omega, \]

where \(\nabla^{4} \equiv \nabla^{2} \nabla^{2}=\Delta\Delta\) is the biharmonic operator and \(f\) is a prescribed source term. To formulate a complete boundary value problem, the biharmonic equation must be complemented by suitable boundary conditions.

Choice of boundary conditions

As we have a fourth order partial differential equation, we are required to supply two boundary conditions. There are two common sets of conditions that people use for the biharmonic equation, namely the clamped condition and the simply supported condition, see for instance Gander and Liu [GL17].

\[\begin{split} \begin{aligned} u =0 &\qquad \frac{\partial u}{\partial n} = 0 # \text{ on } \partial\Omega,\\ u =0 &\qquad \Delta u = 0 \text{ on } \partial \Omega. \end{aligned} \end{split}\]

In this demo we will consider the clamped boundary conditions

\[\begin{split} \begin{aligned} u &= g_D \text{ on } \partial\Omega,\\ \frac{\partial u}{\partial n} &= g_N \text{ on } \partial\Omega,\\ \end{aligned} \end{split}\]

as the simply supported boundary conditions reduces the system into two sequential Poisson problems, named the Ciarlet-Raviart formulation [CR74].

Weak formulation

Multiplying the biharmonic equation by a test function and integrating by parts twice leads to a problem of second-order derivatives, which would require \(H^{2}\) conforming (roughly \(C^{1}\) continuous) basis functions. To solve the biharmonic equation using Lagrange finite element basis functions, the biharmonic equation can be split into two second- order equations (see the Mixed Poisson demo for an example of a mixed method), or a variational formulation can be constructed that imposes weak continuity of normal derivatives between finite element cells. This demo uses a discontinuous Galerkin approach to impose continuity of the normal derivative weakly, see for instance Babuška and Zlámal [BZ73] or Georgoulis and Houston [GH09].

In this demo, we consider equation 3.20-3.21 from Georgoulis and Houston [GH09], a \(\mathcal{C}^0\)-interior penalty method. However, instead of using a broken (discontinuous) finite element space for the unknown \(u\), we use a continuous space, which simplifies the formulation to a weak formulation of the biharmonic problem reads: find \(u \in V_{g_D}\) such that

\[ a(u,v)=L(v) \quad \forall \ v \in V, \]

where the bilinear form is

\[\begin{split} \begin{aligned} a(u, v) &= \sum_{K \in \mathcal{T}} \int_{K} \Delta u \Delta v ~{\rm d}x \\ &+\sum_{E \in \mathcal{E}_h^{\rm int}}\int_{E}\left( - \left<\Delta u \right>[\!\![ \nabla v ]\!\!] - [\!\![ \nabla u ]\!\!] \left<\nabla^{2} v \right> + \frac{\beta}{h_E} [\!\![ \nabla u ]\!\!] [\!\![ \nabla v ]\!\!] \right)~{\rm d}s\\ &+\sum_{E \in \mathcal{E}_h^{\rm ext}}\int_{E}\left( - \Delta u \nabla v \cdot n - \Delta v \nabla u \cdot n + \frac{\beta}{h_E} \nabla u \cdot n \nabla v \cdot n \right)~{\rm d}s \end{aligned} \end{split}\]

and the linear form is

\[ L(v) = \int_{\Omega} fv ~{\rm d}x +\sum_{E \in \mathcal{E}_h^{\rm ext}}\int_{E}\left( -g_N \Delta v + \frac{\beta}{h_E} g_N \nabla v \cdot n \right) ~{\rm d}s. \]

Furthermore, \(\left< u \right> = \frac{1}{2} (u_{+} + u_{-})\), \([\!\![ w ]\!\!] = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}\), \(\beta \ge 0\) is a penalty parameter and \(h_E\) is a measure of the cell size. We use the penalty parameter described in Bringmann et al. [BCS24].

where \(K\) is an element of the mesh, while \(\mathcal{E}_h^{\rm int}\) is the collection of all interior facets, while \(\mathcal{E}_h^{\rm ext}\) is the set of exterior facets.

Note

Note that the Dirichlet condition for \(u\) will be enforced strongly in this example.

Implementation

We follow the example of Georgoulis and Houston [GH09] and use the method of manufactured solutions to construct a \(f\), \(g_D\), \(g_N\) that satisfies

\[ u(x, y) = \sin (2\pi x) \sin (2 \pi y) \text{ in } [0, 1] \times [0, 1]. \]

We first import the modules and functions that the program uses:

from pathlib import Path

from mpi4py import MPI

import numpy as np

import ufl
from dolfinx import default_real_type, default_scalar_type, fem, has_adios2, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, GhostMode

if np.issubdtype(default_real_type, np.float32):  # type: ignore
    print("float32 not yet supported for this demo.")
    exit(0)

We begin by using create_rectangle to create a rectangular Mesh of the domain, and creating a finite element FunctionSpace \(V\) on the mesh.

N = 32
msh = mesh.create_unit_square(
    MPI.COMM_WORLD,
    N,
    N,
    cell_type=CellType.triangle,
    ghost_mode=GhostMode.shared_facet,
)

As noted in Georgoulis and Houston [GH09] second order Lagrange elements yield sub-optimal convergence, and therefore we choose to use third order elements

degree = 3
V = fem.functionspace(msh, ("Lagrange", degree))

The second argument to functionspace is a tuple consisting of (family, degree), where family is the finite element family, and degree specifies the polynomial degree. in this case V consists of third-order, continuous Lagrange finite element functions. For further details of how one can specify finite elements as tuples, see ElementMetaData.

Next, we locate the mesh facets that lie on the boundary \(\Gamma_D = \partial\Omega\). We do this using using exterior_facet_indices which returns all mesh boundary facets (Note: if we are only interested in a subset of those, consider locate_entities_boundary).

tdim = msh.topology.dim
msh.topology.create_connectivity(tdim - 1, tdim)
facets = mesh.exterior_facet_indices(msh.topology)

We now find the degrees-of-freedom that are associated with the boundary facets using locate_dofs_topological

fdim = tdim - 1
dofs = fem.locate_dofs_topological(V=V, entity_dim=fdim, entities=facets)

and use dirichletbc to create a DirichletBC class that represents the boundary condition.

We define the manufactured solution and interpolate it into the function space of our unknown to apply it as a strong boundary condition

x = ufl.SpatialCoordinate(msh)
u_ex = ufl.sin(2 * ufl.pi * x[0]) ** 2 * ufl.sin(2 * ufl.pi * x[1]) ** 2
g_D_expr = fem.Expression(u_ex, V.element.interpolation_points)
g_D = fem.Function(V)
g_D.interpolate(g_D_expr)
bc = fem.dirichletbc(value=g_D, dofs=dofs)

Next, we express the variational problem using UFL.

First, the penalty parameter \(\beta\) is defined. In addition, we define a variable h for the cell diameter \(h_E\), a variable nfor the outward-facing normal vector \(n\) and a variable h_avg for the average size of cells sharing a facet \(\left< h \right> = \frac{1}{2} (h_{+} + h_{-})\). Here, the UFL syntax ('+') and ('-') restricts a function to the ('+') and ('-') sides of a facet.

a = fem.Constant(msh, default_scalar_type(4.0))
vol = ufl.CellVolume(msh)
h = ufl.CellDiameter(msh)
beta = 3.0 * a * degree * (degree - 1) / 8.0 * h("+") ** 2 * ufl.avg(1 / vol)
beta_boundary = 3.0 * a * degree * (degree - 1) * h**2 / vol
n = ufl.FacetNormal(msh)
h_avg = (h("+") + h("-")) / 2.0

After that, we can define the variational problem consisting of the bilinear form \(a\) and the linear form \(L\). We use ufl to derive the manufactured \(f\) and \(g_N\). Note that with dS, integration is carried out over all the interior facets \(\mathcal{E}_h^{\rm int}\), whereas with ds it would be only the facets on the boundary of the domain, i.e. \(\partial\Omega\). The jump operator \([\!\![ w ]\!\!] = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}\) w.r.t. the outward-facing normal vector \(n\) is in UFL available as jump(w, n).

# Define variational problem
u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
f = ufl.div(ufl.grad(ufl.div(ufl.grad(u_ex))))
g_N = ufl.dot(ufl.grad(u_ex), n)

a = (
    ufl.inner(ufl.div(ufl.grad(u)), ufl.div(ufl.grad(v))) * ufl.dx
    - ufl.inner(ufl.avg(ufl.div(ufl.grad(u))), ufl.jump(ufl.grad(v), n)) * ufl.dS
    - ufl.inner(ufl.jump(ufl.grad(u), n), ufl.avg(ufl.div(ufl.grad(v)))) * ufl.dS
    + beta / h_avg * ufl.inner(ufl.jump(ufl.grad(u), n), ufl.jump(ufl.grad(v), n)) * ufl.dS
    - ufl.inner(ufl.div(ufl.grad(u)), ufl.dot(ufl.grad(v), n)) * ufl.ds
    - ufl.inner(ufl.dot(ufl.grad(u), n), ufl.div(ufl.grad(v))) * ufl.ds
    + beta_boundary / h * ufl.inner(ufl.dot(ufl.grad(u), n), ufl.dot(ufl.grad(v), n)) * ufl.ds
)
L = (
    ufl.inner(f, v) * ufl.dx
    - ufl.inner(g_N, ufl.div(ufl.grad(v))) * ufl.ds
    + beta_boundary / h * ufl.inner(g_N, ufl.dot(ufl.grad(v), n)) * ufl.ds
)

We create a LinearProblem object that brings together the variational problem, the Dirichlet boundary condition, and which specifies the linear solver. In this case we use a direct (LU) solver. The solve will compute a solution.

problem = LinearProblem(
    a,
    L,
    bcs=[bc],
    petsc_options_prefix="demo_biharmonic_",
    petsc_options={
        "ksp_type": "preonly",
        "pc_type": "lu",
    },
)
uh = problem.solve()

assert isinstance(uh, fem.Function)
assert problem.solver.getConvergedReason() > 0

We compute the relative \(L^2\)-error between the computed and exact solution:

error = fem.form(ufl.inner(uh - g_D, uh - g_D) * ufl.dx)
local_error = fem.assemble_scalar(error)
glob_error = error.mesh.comm.allreduce(local_error, op=MPI.SUM)
ref_scale = fem.form(ufl.inner(u_ex, u_ex) * ufl.dx)
local_scale = fem.assemble_scalar(ref_scale)
scale = ref_scale.mesh.comm.allreduce(local_scale, op=MPI.SUM)
print("||u_h-u_ex||_{L^2}^2/||u_ex||_L^2^2: ", f"{glob_error / scale:.5e}")
assert glob_error.real / scale.real < 1e-6

The solution can be written to a VTX-file using VTXWriter which can be opened with ParaView

if has_adios2:
    out_folder = Path("out_biharmonic")
    out_folder.mkdir(parents=True, exist_ok=True)
    with io.VTXWriter(msh.comm, out_folder / "biharmonic.bp", [uh]) as file:
        file.write(0.0)

and displayed using pyvista.

try:
    import pyvista

    cells, types, x = plot.vtk_mesh(V)
    grid = pyvista.UnstructuredGrid(cells, types, x)
    grid.point_data["u"] = uh.x.array.real
    grid.set_active_scalars("u")
    plotter = pyvista.Plotter()
    plotter.add_mesh(grid, show_edges=True)
    warped = grid.warp_by_scalar()
    plotter.add_mesh(warped)
    if pyvista.OFF_SCREEN:
        plotter.screenshot(out_folder / "uh_biharmonic.png")
    else:
        plotter.show()
except ModuleNotFoundError:
    print("'pyvista' is required to visualise the solution")
    print("Install 'pyvista' with pip: 'python3 -m pip install pyvista'")

References

[BZ73]

Ivo Babuška and Miloš Zlámal. Nonconforming elements in the finite element method with penalty. SIAM Journal on Numerical Analysis, 10(5):863–875, 1973. URL: http://www.jstor.org/stable/2156320 (visited on 2026-02-25).

[BCS24]

Philipp Bringmann, Carsten Carstensen, and Julian Streitberger. Local parameter selection in the c0 interior penalty method for the biharmonic equation. Journal of Numerical Mathematics, 32(3):257–273, 2024. doi:10.1515/jnma-2023-0028.

[CR74]

P.G. CIARLET and P.A. RAVIART. A mixed finite element method for the biharmonic equation. In Carl de Boor, editor, Mathematical Aspects of Finite Elements in Partial Differential Equations, pages 125–145. Academic Press, 1974. doi:10.1016/B978-0-12-208350-1.50009-1.

[GL17]

Martin J. Gander and Yongxiang Liu. On the definition of dirichlet and neumann conditions for the biharmonic equation and its impact on associated schwarz methods. In Chang-Ock Lee, Xiao-Chuan Cai, David E. Keyes, Hyea Hyun Kim, Axel Klawonn, Eun-Jae Park, and Olof B. Widlund, editors, Domain Decomposition Methods in Science and Engineering XXIII, 303–311. Cham, 2017. Springer International Publishing. doi:10.1007/978-3-319-52389-7_31.

[GH09] (1,2,3,4)

Emmanuil H. Georgoulis and Paul Houston. Discontinuous galerkin methods for the biharmonic problem. IMA Journal of Numerical Analysis, 29(3):573–594, 2009. doi:10.1093/imanum/drn015.