1. Biharmonic equation¶

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 in this demo will look as follows:

1.1. Equation and problem definition¶

The biharmonic equation is a fourth-order elliptic equation. On the domain , , it reads

where is the biharmonic operator and 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 conforming (roughly 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 union of interior facets is denoted by . 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 reads: find such that

where , , is a penalty term and is a measure of the cell size. For the implementation, it is useful to identify the bilinear form

and the linear form

The input parameters for this demos are defined as follows:

• (a unit square)
• (penalty parameter)
• (source term)

1.2. 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.

1.2.1. UFL form file¶

First we define the variational problem in UFL in the file called Biharmonic.ufl.

In the UFL file, the finite element space is defined:

# Elements
element = FiniteElement("Lagrange", triangle, 2)


On the space element, trial and test functions, and the source term are defined:

# Trial and test functions
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)


Next, the outward unit normal to cell boundaries and a measure of the cell size are defined. The average size of cells sharing a facet will be used (h_avg). The UFL syntax ('+') and ('-') restricts a function to the ('+') and ('-') sides of a facet, respectively. The penalty parameter alpha is made a Constant so that it can be changed in the program without regenerating the code.

# Normal component, mesh size and right-hand side
n  = element.cell().n
h_avg = (h('+') + h('-'))/2

# Parameters
alpha = Constant(triangle)


Finally the bilinear and linear forms are defined. Integrals over internal facets are indicated by *dS.

# Bilinear form

# Linear form
L = f*v*dx


1.2.2. 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 , 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()
{
// Create mesh
UnitSquare mesh(32, 32);


The source function, a function for the cell size and the penalty term are declared:

// Create functions
Source f;
Constant alpha(8.0);


A function space object, which is defined in the generated code, is created:

// Create function space
Biharmonic::FunctionSpace V(mesh);


The Dirichlet boundary condition on 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
Constant u0(0.0);
DirichletBoundary boundary;
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 plotted to the screen and written to a file in VTK format:

  // Plot solution
plot(u);

// Save solution in VTK format
File file("biharmonic.pvd");
file << u;

return 0;
}


1.3. Complete code¶

1.3.1. Complete UFL file¶

# Elements
element = FiniteElement("Lagrange", triangle, 2)

# Trial and test functions
u = TrialFunction(element)
v = TestFunction(element)
f = Coefficient(element)

# Normal component, mesh size and right-hand side
n  = element.cell().n
h_avg = (h('+') + h('-'))/2

# Parameters
alpha = Constant(triangle)

# Bilinear form

# Linear form
L = f*v*dx


1.3.2. Complete main file¶

#include <dolfin.h>
#include "Biharmonic.h"

using namespace dolfin;

// 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]);
}

};

// Sub domain for Dirichlet boundary condition
class DirichletBoundary : public SubDomain
{
bool inside(const Array<double>& x, bool on_boundary) const
{
return on_boundary;
}
};

int main()
{
// Create mesh
UnitSquare mesh(32, 32);

// Create functions
Source f;
Constant alpha(8.0);

// Create function space
Biharmonic::FunctionSpace V(mesh);

// Define boundary condition
Constant u0(0.0);
DirichletBoundary boundary;
DirichletBC bc(V, u0, boundary);

// Define variational problem
Biharmonic::BilinearForm a(V, V);
Biharmonic::LinearForm L(V);
a.alpha = alpha; L.f = f;

// Compute solution
Function u(V);
solve(a == L, u, bc);

// Plot solution
plot(u);

// Save solution in VTK format
File file("biharmonic.pvd");
file << u;

return 0;
}