demo_biharmonic.py.rst

Biharmonic equation

This demo is implemented in a single Python file, demo_biharmonic.py, which contains both the variational forms and the solver.

This demo illustrates how to:

• Solve a linear partial differential equation
• Use a discontinuous Galerkin method
• Solve a fourth-order differential equation

The solution for u in this demo will look as follows:

Equation and problem definition

The biharmonic equation is a fourth-order elliptic equation. On the domain Ω ⊂ ℝ^d, 1 ≤ d ≤ 3, it reads

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

where ∇⁴ ≡ ∇²∇² 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.

Multiplying the biharmonic equation by a test function and integrating by parts twice leads to a problem second-order derivatives, which would requires H² conforming (roughly C¹ 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. The demo uses a discontinuous Galerkin approach to impose continuity of the normal derivative weakly.

Consider a triangulation 𝒯 of the domain Ω, 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

V = {v∈H₀¹(Ω) : v∈P_k(K) ∀ K∈𝒯}

and considering the boundary conditions

$$\begin{aligned} u &= 0 \quad {\rm on} \ \partial\Omega \\\ \nabla^{2} u &= 0 \quad {\rm on} \ \partial\Omega \end{aligned}$$

a weak formulation of the biharmonic problem reads: find u ∈ V such that

a(u, v) = L(v) ∀ v ∈ V,

where the bilinear form is

$$a(u, v) = \sum_{K \in \mathcal{T}} \int_{K} \nabla^{2} u \nabla^{2} v \, {\rm d}x \\ +\sum_{E \in \mathcal{E}_h^{\rm int}}\left(\int_{E} \frac{\alpha}{h_E} [\!\![ \nabla u ]\!\!] [\!\![ \nabla v ]\!\!] \, {\rm d}s - \int_{E} \left<\nabla^{2} u \right>[\!\![ \nabla v ]\!\!] \, {\rm d}s - \int_{E} [\!\![ \nabla u ]\!\!] \left<\nabla^{2} v \right> \, {\rm d}s\right)$$

and the linear form is

$$L(v) = \int_{\Omega} fv \, {\rm d}x$$

Furthermore, $\left&lt; u \right&gt; = \frac{1}{2} (u_{+} + u_{-})$, [​​[w]​​] = w₊ ⋅ n₊ + w₋ ⋅ n₋, α ≥ 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:

• Ω = [0, 1] × [0, 1] (a unit square)
• α = 8.0 (penalty parameter)
• f = 4.0π⁴sin (πx)sin (πy) (source term)

Implementation

This demo is implemented in the demo_biharmonic.py file.

First, the :pydolfin module is imported:

from dolfin import *

Specifying matplotlib as backend:

%matplotlib inline
import matplotlib.pyplot as plt
parameters["plotting_backend"]="matplotlib"

Next, some parameters for the form compiler are set:

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["optimize"] = True

A mesh is created, and a quadratic finite element function space:

# Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet"

# Create mesh and define function space
mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 2)

A subclass of :pySubDomain <dolfin.cpp.SubDomain>, DirichletBoundary is created for later defining the boundary of the domain:

# Define Dirichlet boundary
class DirichletBoundary(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary

A subclass of :pyExpression <dolfin.functions.expression.Expression>, Source is created for the source term f:

class Source(Expression):
    def eval(self, values, x):
        values[0] = 4.0*pi**4*sin(pi*x[0])*sin(pi*x[1])

The Dirichlet boundary condition is created:

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, DirichletBoundary())

On the finite element space V, trial and test functions are created:

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)

A function for the cell size h is created, as is a function for the average size of cells that share a facet (h_avg). The UFL syntax ('+') and ('-') restricts a function to the ('+') and ('-') sides of a facet, respectively. The unit outward normal to cell boundaries (n) is created, as is the source term f and the penalty parameter alpha. The penalty parameters is made a :pyConstant <dolfin.functions.constant.Constant> so that it can be changed without needing to regenerate code. :

# Define normal component, mesh size and right-hand side
h = CellSize(mesh)
h_avg = (h('+') + h('-'))/2.0
n = FacetNormal(mesh)
f = Source(degree=2)

# Penalty parameter
alpha = Constant(8.0)

The bilinear and linear forms are defined:

# Define bilinear form
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

# Define linear form
L = f*v*dx

A :pyFunction <dolfin.functions.function.Function> is created to store the solution and the variational problem is solved:

# Solve variational problem
u = Function(V)
solve(a == L, u, bc)

The computed solution is written to a file in VTK format and plotted to the screen. :

# Save solution to file
file = File("biharmonic.pvd")
file << u

# Plot solution
plt.figure()
plot(u, interactive=True)