Biharmonic equation (C++)¶
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 \(\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.
Multiplying the biharmonic equation by a test function and integrating by parts twice leads to a problem second-order derivatives, which would requires \(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. The 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¶
The implementation is split in two files, a form file containing the definition of the variational forms expressed in UFL and the solver which is implemented in a C++ file.
Running this demo requires the files: main.cpp
,
Biharmonic.ufl
and CMakeLists.txt
.
UFL form file¶
The UFL file is implemented in Biharmonic.ufl
, and the
explanation of the UFL file can be found at here.
C++ program¶
The DOLFIN interface and the code generated from the UFL input is included, and the DOLFIN namespace is used:
#include <dolfin.h>
#include "Biharmonic.h"
using namespace dolfin;
A class Source
is defined for the function \(f\), with the
function Expression::eval
overloaded:
// Source term
class Source : public Expression
{
public:
void eval(Array<double>& values, const Array<double>& x) const
{
values[0] = 4.0*std::pow(DOLFIN_PI, 4)*
std::sin(DOLFIN_PI*x[0])*std::sin(DOLFIN_PI*x[1]);
}
};
A boundary subdomain is defined, which in this case is the entire boundary:
// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{ return on_boundary; }
};
The main part of the program is begun, and a mesh is created with 32 vertices in each direction:
int main()
{
// Make mesh ghosted for evaluation of DG terms
parameters["ghost_mode"] = "shared_facet";
// Create mesh
auto mesh = std::make_shared<UnitSquareMesh>(32, 32);
The source function, a function for the cell size and the penalty term are declared:
// Create functions
auto f = std::make_shared<Source>();
auto alpha = std::make_shared<Constant>(8.0);
A function space object, which is defined in the generated code, is created:
// Create function space
auto V = std::make_shared<Biharmonic::FunctionSpace>(mesh);
The Dirichlet boundary condition on \(u\) is constructed by
defining a Constant
which is equal to zero, defining the
boundary (DirichletBoundary
), and using these, together with
V
, to create bc
:
// Define boundary condition
auto u0 = std::make_shared<Constant>(0.0);
auto boundary = std::make_shared<DirichletBoundary>();
DirichletBC bc(V, u0, boundary);
Using the function space V
, the bilinear and linear forms are
created, and function are attached:
// Define variational problem
Biharmonic::BilinearForm a(V, V);
Biharmonic::LinearForm L(V);
a.alpha = alpha; L.f = f;
A Function
is created to hold the solution and the
problem is solved:
// Compute solution
Function u(V);
solve(a == L, u, bc);
The solution is then written to a file in VTK format and plotted to the screen:
// Save solution in VTK format
File file("biharmonic.pvd");
file << u;
// Plot solution
plot(u);
interactive();
return 0;
}