Biharmonic equation
This demo is implemented in a single Python file,
demo_biharmonic.py
, which contains both the variational forms
and the solver. It illustrates how to:
Solve a linear partial differential equation
Use a discontinuous Galerkin method
Solve a fourth-order differential equation
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
where \(\nabla^{4} \equiv \nabla^{2} \nabla^{2}\) 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.
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 a mixed method for the Poisson equation), 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.
Consider a triangulation \(\mathcal{T}\) of the domain \(\Omega\), where the set of interior facets is denoted by \(\mathcal{E}_h^{\rm int}\). Functions evaluated on opposite sides of a facet are indicated by the subscripts \(+\) and \(-\). Using the standard continuous Lagrange finite element space
and considering the boundary conditions
a weak formulation of the biharmonic problem reads: find \(u \in V\) such that
where the bilinear form is
and the linear form is
Furthermore, \(\left< u \right> = \frac{1}{2} (u_{+} + u_{-})\), \([\!\![ w ]\!\!] = w_{+} \cdot n_{+} + w_{-} \cdot n_{-}\), \(\alpha \ge 0\) is a penalty parameter and \(h_E\) is a measure of the cell size.
The input parameters for this demo are defined as follows:
\(\Omega = [0,1] \times [0,1]\) (a unit square)
\(\alpha = 8.0\) (penalty parameter)
\(f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)\) (source term)
Implementation
We first import the modules and functions that the program uses:
from mpi4py import MPI
from petsc4py.PETSc import ScalarType # type: ignore
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, GhostMode
We begin by using create_rectangle
to create a rectangular
Mesh
of the domain, and creating a
finite element FunctionSpace
\(V\) on the mesh.
msh = mesh.create_rectangle(
comm=MPI.COMM_WORLD,
points=((0.0, 0.0), (1.0, 1.0)),
n=(32, 32),
cell_type=CellType.triangle,
ghost_mode=GhostMode.shared_facet,
)
V = fem.functionspace(msh, ("Lagrange", 2))
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
second-order, continuous Lagrange finite element functions.
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
dofs = fem.locate_dofs_topological(V=V, entity_dim=1, entities=facets)
and use dirichletbc
to create a
DirichletBC
class that represents the boundary condition. In this case, we impose
Dirichlet boundary conditions with value \(0\) on the entire boundary
\(\partial\Omega\).
bc = fem.dirichletbc(value=ScalarType(0), dofs=dofs, V=V)
Next, we express the variational problem using UFL.
First, the penalty parameter \(\alpha\) is defined. In addition, we define
a variable h
for the cell diameter \(h_E\), a variable n
for 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.
alpha = ScalarType(8.0)
h = ufl.CellDiameter(msh)
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\). The source term is prescribed
as \(f = 4.0 \pi^4\sin(\pi x)\sin(\pi y)\). 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)
x = ufl.SpatialCoordinate(msh)
f = 4.0 * ufl.pi**4 * ufl.sin(ufl.pi * x[0]) * ufl.sin(ufl.pi * x[1])
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
+ alpha / h_avg * ufl.inner(ufl.jump(ufl.grad(u), n), ufl.jump(ufl.grad(v), n)) * ufl.dS
)
L = ufl.inner(f, v) * ufl.dx
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
The solution can be written to a XDMFFile
file visualization with ParaView or VisIt
with io.XDMFFile(msh.comm, "out_biharmonic/biharmonic.xdmf", "w") as file:
V1 = fem.functionspace(msh, ("Lagrange", 1))
u1 = fem.Function(V1)
u1.interpolate(uh)
file.write_mesh(msh)
file.write_function(u1)
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:
pyvista.start_xvfb(wait=0.1)
plotter.screenshot("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'")