1. Introduction
  2. Variational Formulation
  3. Commented Program
  4. Results
  5. Complete Source Code
@section intro Introduction
In this tutorial the implementation of inhomogeneous Dirichlet boundary 
conditions using a Dirichlet lift ansatz is shown.

The equation solved is the reaction-diffusion equation
\f[ - \Delta u + u = f \f]
in the unit square \f$ \Omega = (0,1)^2 \f$ with source term
f(x,y) = -5\exp\left(-(x-1/2)^2-(y-1/2)^2\right)
and Dirichlet boundary conditions
u = 
\sin \pi x, &(x,y)\in\Gamma_{\rm inh}=\{(x,y)\in\partial\Omega\mid y=0\},\\
0,          &(x,y)\in\partial\Omega\setminus\Gamma_{\rm inh}.

A suitable Dirichlet lift is
\f[ g = \sin \pi x \; \cos \frac{\pi}{2} y. \f]
By setting \f$ u = \tilde{u} + g \f$ the reaction-diffusion equation simplifies
to a problem in \f$ \tilde{u} \f$
- \Delta \tilde{u} + \tilde{u} = f + \Delta g - g, \qquad x\in\Omega,
with homogeneous Dirichlet boundary conditions
\tilde{u} = 0, \qquad x\in\partial\Omega.

@section variation Variational Formulation
Now we derive the corresponding variational formulation of the introduced problem:
find \f$ \tilde{u}\in H^1_0(\Omega) \f$ such that
\int_{\Omega}\nabla\tilde{u}\cdot\nabla v\;{\rm d}x
+ \int_{\Omega}\tilde{u}v\;{\rm d}x
\int_{\Omega}fv\;{\rm d}x
-\int_{\Omega}\nabla g\cdot\nabla v\;{\rm d}x
- \int_{\Omega}gv\;{\rm d}x \qquad\forall v\in H^1_0(\Omega).

@section commented Commented Program
First, system files
@until iostream
and concepts files
@until toolbox
are included. With the following using directive
@until using
we do not need to prepend <tt>concepts::</tt> to <tt>Real</tt> everytime.

@subsection main Main Program
The mesh is read from three files containing the coordinates, the
elements and the boundary attributes of the mesh.
@until endl
<!--Please see the tutorial @ref mesh "Import 2D Mesh" for a concise explanation of how to write these import files and define the attributes for edges etc. -->
In our example the edges are given the following attributes: 1 (bottom), 2 (right), 3 (top) and 4 (left).

Now the mesh is plotted using a scaling factor of 100, a greyscale of 1.0 
and one point per edge.
@until drawMeshEPS

The Dirichlet lift and its gradient are defined.
@skip concepts
@until DirichletLiftGrady

The problem for \f$ \tilde{u} \f$ has homogeneous Dirichlet boundary conditions at 
all four edges.
@skip concepts
@until bc.add

Using the mesh and the homogeneous Dirichlet boundary conditions, the space 
can be built. We refine the space two times and set the polynomial degree to three.
Then the elements of the space are built and the space is plotted.
@skip hp2D
@until drawMeshEPS

The right hand side is computed, containing the integrals of the source term, 
the Dirichlet lift and its gradient.
@skip hp2D
@until endl

The system matrix is computed from the bilinear form \c hp2D::Laplace<Real> and 
the bilinear form \c hp2D::Identity<Real>.
@skip hp2D
@until endl

We solve the equation using the conjugate gradient method with a tolerance of 
1e-6 and a maximum number of iterations of 200.
@skip concepts
@until endl

In order to add the Dirichlet lift to the solution of the homogeneous Dirichlet 
problem we transform the solution, that is given as a vector related to the basis
functions of \c spc, into an element formula, that can be evaluated at each point
in every cell.
@skip concepts
@until endl

We plot the solution, the homogeneous Dirichlet solution and the Dirichlet lift
using <tt>graphics::MatlabGraphics</tt>. To this end the shape functions are 
computed on equidistant points using the trapezoidal quadrature rule.
@skip hp2D
@until MatlabGraphics
@until MatlabGraphics
@until MatlabGraphics

Finally, exceptions are caught and a sane return value is given back.
@until }
@until }

@section results Results
Output of the program:

Mesh: Import2dMesh(ncell = 1)

RHS Vector: Vector(121, [-7.813458e-01, -1.331759e-01, -1.446913e-04, -8.977796e-02, -1.975711e-02, -1.504745e-02, -1.103641e-04, -3.490666e-03, 4.167881e-05, -1.009906e+00, -1.735826e-01, 3.466376e-04, -1.628611e-01, 8.126026e-03, -2.796551e-02, 4.594644e-05, -1.436936e-03, 1.759809e-05, -8.586824e-01, -6.786590e-01, -1.598694e-01, 5.110721e-03, -1.388192e-01, -6.395889e-03, -1.248594e-01, -3.466673e-03, -2.580936e-02, 8.124231e-04, -1.244348e-03, 5.860425e-05, -8.122847e-02, -1.558898e-02, -1.461248e-02, 2.876690e-04, -3.028574e-03, 1.411786e-04, -7.813458e-01, -1.331759e-01, -1.446913e-04, -1.628611e-01, 8.126026e-03, -2.796551e-02, 4.594644e-05, 1.436936e-03, -1.759809e-05, -8.977796e-02, -1.975711e-02, -1.504745e-02, -1.103641e-04, 3.490666e-03, -4.167881e-05, -6.786590e-01, -8.122847e-02, -1.558898e-02, -1.248594e-01, -3.466673e-03, -1.461248e-02, 2.876690e-04, 3.028574e-03, -1.411786e-04, -1.388192e-01, -6.395889e-03, -2.580936e-02, 8.124231e-04, 1.244348e-03, -5.860425e-05, -4.822921e-01, -5.869802e-01, -9.922001e-02, 6.668103e-03, -9.533732e-02, 3.713134e-03, -1.236099e-01, -9.222017e-03, -2.002021e-02, 1.475922e-03, 8.660092e-04, -9.094635e-05, -6.180870e-02, -9.103448e-03, -1.221454e-02, 6.601945e-04, 2.115183e-03, -2.198685e-04, -6.016121e-02, -8.736539e-03, -8.218690e-03, 9.139240e-04, 8.895484e-04, -2.644021e-04, -7.032423e-02, -1.180407e-02, -1.147942e-02, 1.894019e-03, 3.595180e-04, -1.091848e-04, -4.822921e-01, -9.922001e-02, 6.668103e-03, -6.180870e-02, -9.103448e-03, -1.221454e-02, 6.601945e-04, -2.115183e-03, 2.198685e-04, -9.533732e-02, 3.713134e-03, -2.002021e-02, 1.475922e-03, -8.660092e-04, 9.094635e-05, -6.016121e-02, -8.736539e-03, -1.147942e-02, 1.894019e-03, -3.595180e-04, 1.091848e-04, -8.218690e-03, 9.139240e-04, -8.895484e-04, 2.644021e-04])

System Matrix: SparseMatrix(121x121, HashedSparseMatrix: 1521 (1.038863e+01%) entries bound.)

Solver: CG(solves SparseMatrix(121x121, HashedSparseMatrix: 1521 (1.038863e+01%) entries bound.), eps = 6.210533e-13, it = 37, relres = 0)

Solution: FrmE_Sum(ElementFormulaVector<1>(hp2D::Value<double>(), Vector(121, [-5.387654e-01, -3.006218e-01, 1.564124e-02, -1.425659e-01, -4.809499e-03, -1.088256e-01, 1.611971e-02, 8.625588e-03, -6.891406e-03, -7.373583e-01, -4.048944e-01, 1.925875e-02, -1.941122e-01, 2.767690e-03, -1.016416e-01, 3.259733e-03, -1.569518e-03, 2.857100e-04, -8.461932e-01, -6.183710e-01, -2.393483e-01, 9.380839e-03, -2.230022e-01, -2.981412e-03, -1.722831e-01, -6.882467e-03, -6.528885e-02, 2.490903e-03, -1.108972e-03, -1.050398e-05, -1.618117e-01, -6.560393e-03, -3.530271e-02, 1.727733e-03, -4.305751e-03, 3.179923e-04, -5.387655e-01, -3.006218e-01, 1.564124e-02, -1.941122e-01, 2.767687e-03, -1.016416e-01, 3.259733e-03, 1.569518e-03, -2.857100e-04, -1.425659e-01, -4.809502e-03, -1.088256e-01, 1.611971e-02, -8.625588e-03, 6.891407e-03, -6.183710e-01, -1.618117e-01, -6.560394e-03, -1.722831e-01, -6.882469e-03, -3.530271e-02, 1.727733e-03, 4.305752e-03, -3.179924e-04, -2.230022e-01, -2.981412e-03, -6.528885e-02, 2.490903e-03, 1.108972e-03, 1.050400e-05, -4.161721e-01, -5.639852e-01, -1.159405e-01, 2.828855e-03, -1.454700e-01, 1.444007e-03, -1.596678e-01, -4.353042e-03, -4.288914e-02, 1.375478e-03, 4.954845e-04, -1.172436e-04, -1.224176e-01, -1.613845e-03, -2.602445e-02, -1.262254e-04, 2.824662e-03, 6.029993e-05, -1.041387e-01, 1.097789e-03, -7.650684e-02, -1.372658e-02, -1.375732e-02, -7.269106e-03, -1.270253e-01, -1.308806e-03, -2.361722e-02, 2.517779e-03, -5.561139e-04, 1.292617e-04, -4.161721e-01, -1.159405e-01, 2.828853e-03, -1.224176e-01, -1.613847e-03, -2.602445e-02, -1.262255e-04, -2.824662e-03, -6.029997e-05, -1.454700e-01, 1.444005e-03, -4.288914e-02, 1.375478e-03, -4.954845e-04, 1.172437e-04, -1.041387e-01, 1.097788e-03, -2.361722e-02, 2.517779e-03, 5.561139e-04, -1.292618e-04, -7.650684e-02, -1.372658e-02, 1.375732e-02, 7.269106e-03])) + ParsedFormula<Real>((sin(pi*x)*cos(pi/2*y))))

Matlab plot of the solution: Matlab output

@section complete Complete Source Code
@author Dirk Klindworth, 2011
#include <iostream>
#include "basics.hh"
#include "formula.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "hp2D.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"
int main(int argc, char** argv) {
try {
concepts::Import2dMesh msh(SOURCEDIR "/applications/unitsquareCoordinates.dat",
SOURCEDIR "/applications/unitsquareElements.dat",
SOURCEDIR "/applications/unitsquareEdges.dat");
std::cout << std::endl << "Mesh: " << msh << std::endl;
graphics::drawMeshEPS(msh, "mesh.eps",100,1.0,1);
// ** define Dirichlet lift **
concepts::ParsedFormula<> DirichletLift("(sin(pi*x)*cos(pi/2*y))");
concepts::ParsedFormula<> DirichletLiftGradx("(pi*cos(pi*x)*cos(pi/2*y))");
concepts::ParsedFormula<> DirichletLiftGrady("(-pi/2*sin(pi*x)*sin(pi/2*y))");
// ** homogeneous Dirichlet boundary conditions **
for (uint i = 1; i <= 4; ++i)
// ** space **
hp2D::hpAdaptiveSpaceH1 spc(msh, 2, 3, &bc);
graphics::drawMeshEPS(spc, "space.eps",100,1.0,1);
// ** right hand side with added Dirichlet lift **
hp2D::Riesz<Real> lform(concepts::ParsedFormula<>("(-5*exp(-(x-0.5)*(x-0.5)-(y-0.5)*(y-0.5)))"));
concepts::Vector<Real> rhs(spc, lform);
hp2D::GradLinearForm<Real> lformDirichletGrad(DirichletLiftGradx,DirichletLiftGrady);
concepts::Vector<Real> rhsDirichletGrad(spc, lformDirichletGrad);
hp2D::Riesz<Real> lformDirichlet(DirichletLift);
concepts::Vector<Real> rhsDirichlet(spc, lformDirichlet);
rhs = rhs - rhsDirichletGrad - rhsDirichlet;
std::cout << std::endl << "RHS Vector:" << std::endl << rhs << std::endl;
// ** left hand side matrix **
M.addInto(A, 1.0);
std::cout << std::endl << "System Matrix:" << std::endl << A << std::endl;
// ** solve **
concepts::CG<Real> solver(A, 1e-6, 200);
solver(rhs, sol);
std::cout << std::endl << "Solver:" << std::endl << solver << std::endl;
// ** transform solution to formula (in order to add Dirichlet lift) **
std::cout << std::endl << "Solution:" << std::endl << solFormula + DirichletLift << std::endl;
// ** plot solution **
hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, 8);
graphics::MatlabGraphics(spc, "inhomDirichlet", solFormula + DirichletLift);
graphics::MatlabGraphics(spc, "inhomDirichlet0", solFormula);
graphics::MatlabGraphics(spc, "inhomDirichletLift", DirichletLift);
std::cout << e << std::endl;
return 1;
return 0;
Draws a picture of data in Matlab format and stores the result in a single file.
Definition: matlab.hh:114
The approximated function in a FE space.
Definition: function.hh:33
Solves a symmetric system of linear equations with conjugate gradients (CG).
Definition: cg.hh:39
static std::unique_ptr< concepts::QuadRuleFactoryTensor2d > & factory()
Access to the quadrature rule, which is valid for all elements of this type (hp2D::IntegrableQuad).
Definition: quad.hh:145
Definition: defines.hh:13
Base class for exceptions.
Definition: exceptions.hh:86
MeshEPS< Real > drawMeshEPS(concepts::Mesh &msh, std::string filename, const Real scale=100, const Real greyscale=1.0, const unsigned int nPoints=2)
Trampoline function to create a MeshEPS.
void add(const Set< Attribute > &attrib, const Boundary &bcObject)
Adds a boundary condition with this attribute to the list of boundary conditions.
void compress(Real threshold=EPS)
Compresses the matrix by dropping small entries.
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
Vectorial formula created from a FE function.
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.
Linear form in 2D.
Definition: linearForm.hh:138
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:
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich