On this page you find an easy example of an elliptic partial differential equation (PDE) solved with Concepts. Starting with this PDE, we derive its associated variational formulation and solve it numerically. This example is the standard setting in the Concepts tutorials where you can find many extensions to this simple code.
We want to solve the following elliptic PDE with homogeneous Dirichlet boundary conditions
with the right hand side and . This problem has the unique solution .
To solve this problem with Concepts, we first need its variational formulation. Let be the Sobolev space of square integrable and weak differentiable functions with vanishing trace on . Then we seek the weak solution that satisfies
Using integration by parts (Gauss theorem), this is equivalent to
where denotes the unit normal vector outward to . Using that has vanishing trace, we get
This is the variational formulation of our problem. Let now be a finite dimensional subspace of . The Galerkin approximation of our solution in the space is given as the (unique) solution that fulfills the equation
Let be a basis of the space , then we can find a unique vector such that . Plugging this into the equation above and using the linearity of the integral, we obtain the linear equation
For , and and matrices , and a vector this is equivalent to the linear equation
with the stiffness matrix , the mass matrix and the load vector of the right hand side.
Concepts uses the hp-finite element method to approximate the solution of a PDE, which means it calculates the Galerkin approximation of the solution in the finite dimensional subspace
where is a mesh on with mesh size and is a cell in the mesh.
Now let us see how to calculate the FEM-approximation of the solution of the PDE with Concepts. First we have to include some headers to use the functions of specific Concepts packages.
@skip hp2D.hh @until concepts.hh
hp2D includes all functions that are needed to build a representation of the hp-finite element space in two dimensions. The meta-package
concepts includes the data structures for Concepts matrices and linear solvers, allows you to represent the right hand side via its symbolic formula and is used for the graphical output of the mesh and the solution.
allows us to use the datatype
concepts::Real just by typing
Next we start with the actual code, consisting of a main program.
First we build the geometry, a square with one vertex in the origin, and another one in the point . The third parameter is the attribute of all edges, which we set to to .
Now we prescribe homogeneous Dirichlet boundary conditions to all edges of attribute , i.e. to all edges.
With the mesh and the boundary conditions at hand we can build the hp-finite element space.
This is a finite element space on the square, with a twice refined mesh (a single refinement of a two dimensional mesh divides an element in the mesh into four elements, so in our case we have 16 subelements in our square) and polynomial degree six on every element.
Now we define the linear form on the right hand side of the variational formulation.
In the first line we define a representation of and in the second line we define a representation of the linear form .
In the next line we calculate the vector of the right hand side using the space
space and the linear form
We compute the stiffness and mass matrices in a similar way. We first initialize the bilinear forms as
The first one represents the bilinear form that is the variational representation of the Laplace operator and the second
concepts::BilinearForm represents the mapping .
Now we build the matrices using a
concepts::SparseMatrix constructor that computes the matrix belonging the space
space and the bilinear forms
Before we solve the linear equation, we add the mass matrix into the stiffness matrix.
We solve the linear equation using a conjugate gradient (CG) solver which is appropriate as the the system matrix is symmetric. We set the error tolerance to 1e-12 and the maximum number of iterations to 400.
Now we want to save the information about the mesh and the values of the solution on the integration points in a MAT-file that can be read by Matlab or gnu Octave.
First we set the integration rule to the trapezoidal rule and recompute the shape functions on this points.
Then we save the information about the mesh in a MAT-file named "firstSolution.mat" and add the values of the solution on the integration points as a vector named "u".
Finally we just leave the main program.
If you compile the program and run it, you compute a FE-approximation of the solution of the PDE. In the directory where you run the program you will find a MAT-file named "firstSolution.mat". If you use Matlab you can load this MAT-file and plot the solution entering the following instructions in your Matlab command line
This gives you the following visualization of the solution .
You can visualize the vertexes and the edges of the mesh entering the following instructions in your Matlab command line