# howToGetStarted.cc

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.

# Introduction

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.

# Commented Program

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


The package 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.

The instruction

allows us to use the datatype concepts::Real just by typing Real.

Next we start with the actual code, consisting of a main program.

int main(int argc, char** argv) {

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 .

concepts::Square mesh(1, 1, 1);

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.

uint levelOfRefinement = 2;
uint polynomialDegree = 6;
space.rebuild();

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.

concepts::ParsedFormula<> f("((2*pi*pi+1)*sin(pi*x)*sin(pi*y))");

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

concepts::Vector<Real> rhs(space, lform);

We compute the stiffness and mass matrices in a similar way. We first initialize the bilinear forms as concepts::BilinearForm

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 la and id, respectively.

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.

std::unique_ptr<concepts::Operator<Real> > solver(nullptr);
solver.reset(new concepts::CG<Real>(S, 1e-12, 400));
(*solver)(rhs, sol);

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.

space.recomputeShapefunctions();

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

graphics::MatlabBinaryGraphics mbg(space, "firstSolution");

Finally we just leave the main program.

return 0;
}

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

> patch(x(msh), y(msh), u(msh), 'edgecolor', 'none')
> colormap jet

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

> hold on ; plot(x(vtxmsh), y(vtxmsh), 'ko') ; hold off
> hold on ; plot(x(edgemsh), y(edgemsh), 'k-') ; hold off

# Complete Source Code

#include "hp2D.hh"
#include "concepts.hh"
int main(int argc, char** argv) {
// build a mesh of a unit square with edges of attribute 1
concepts::Square mesh(1, 1, 1);
// set homogeneous Dirichlet boundary conditions for every edge of attribute 1 (i.e. all edges)
// build the finite element space V_2_5
uint levelOfRefinement = 2;
uint polynomialDegree = 6;
space.rebuild();
// define the linear form f of the right hand side
concepts::ParsedFormula<> f("((2*pi*pi+1)*sin(pi*x)*sin(pi*y))");
// compute the vector for the finite element discretization
concepts::Vector<Real> rhs(space, lform);
// define the bilinear forms
// compute the mass and stiffness matrices
// add mass matrix to stiffness matrix
// compute the solution of the linear problem using the CG-method
std::unique_ptr<concepts::Operator<Real> > solver(nullptr);
solver.reset(new concepts::CG<Real>(S, 1e-12, 400));
(*solver)(rhs, sol);
// set integration points to be equidistant using trapezoidal rule
// save graphical output of the solution in a MAT-file
graphics::MatlabBinaryGraphics mbg(space, "firstSolution");
return 0;
}
Solves a symmetric system of linear equations with conjugate gradients (CG).
Definition: cg.hh:39
static std::unique_ptr< concepts::QuadRuleFactoryTensor2d > & factory()
@ TRAPEZE
Definition: defines.hh:13
void add(const Set< Attribute > &attrib, const Boundary &bcObject)
Adds a boundary condition with this attribute to the list of boundary conditions.
void addInto(Matrix< H > &dest, const I fact, const uint rowoffset=0, const uint coloffset=0) const
This matrix is added as block to the given matrix dest.
A function class to calculate element matrices for the mass matrix.
Definition: bf_identity.hh:58
void addSolution(const concepts::Space< G > &spc, const std::string name, const concepts::Vector< F > &sol, const concepts::ElementFunction< F, G > *fun=0)
Adds a solution vector to the current matfile.
Class that allows to store graphical infomations in .mat files to use them in Matlab.
Class to describe an element of the boundary.
Definition: boundary.hh:35
void rebuild(bool sameIndices=false)
Rebuilds the mesh and the elements due to adjustment orders.
Definition: square.hh:22
Linear form in 2D.
Definition: linearForm.hh:62
A function class to calculate element matrices for the Laplacian.
Definition: bf_laplace.hh:91
virtual void recomputeShapefunctions()
Recompute shape functions, e.g.
Attributes for elements of the topology.
Definition: connector.hh:22
double Real
Type normally used for a floating point number.
Definition: typedefs.hh:36
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich