# 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 numpy as np
import ufl
from dolfinx import fem, io, mesh, plot
from dolfinx.fem.petsc import LinearProblem
from dolfinx.mesh import CellType, GhostMode
from ufl import CellDiameter, FacetNormal, avg, div, dS, dx, grad, inner, jump, pi, sin
```

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 `locate_entities_boundary`

and providing a marker
function that returns `True`

for points `x`

on the boundary and
`False`

otherwise.

```
facets = mesh.locate_entities_boundary(
msh,
dim=1,
marker=lambda x: np.logical_or.reduce(
(np.isclose(x[0], 0.0), np.isclose(x[0], 1.0), np.isclose(x[1], 0.0), np.isclose(x[1], 1.0))
),
)
```

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 = CellDiameter(msh)
n = 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 * pi**4 * sin(pi * x[0]) * sin(pi * x[1])
a = (
inner(div(grad(u)), div(grad(v))) * dx
- inner(avg(div(grad(u))), jump(grad(v), n)) * dS
- inner(jump(grad(u), n), avg(div(grad(v)))) * dS
+ alpha / h_avg * inner(jump(grad(u), n), jump(grad(v), n)) * dS
)
L = inner(f, v) * 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={"ksp_type": "preonly", "pc_type": "lu"})
uh = problem.solve()
```

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'")
```