# Introduction

In this tutorial the implementation of Robin boundary conditions using trace spaces is shown.

The equation solved is the reaction-diffusion equation

in the unit square with inhomogeneous Neumann boundary condition

and Robin boundary condition

where , and .

# Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find such that

# Commented Program

First, system files

#include <iostream>

and concepts files

#include "basics.hh"
#include "formula.hh"
#include "function.hh"
#include "geometry.hh"
#include "graphics.hh"
#include "hp1D.hh"
#include "hp2D.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"

are included. With the following using directive

we do not need to prepend concepts:: to Real everytime.

## Main Program

The mesh is read from three files containing the coordinates, the elements and the boundary attributes of the mesh.

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;

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.

graphics::drawMeshEPS(msh, "mesh.eps",100,1.0,1);

Since we do not have any homogeneous Dirichlet boundary conditions we do not need to set them using concepts::BoundaryConditions bc.

The formula of the inhomogeneous Neumann data is defined

concepts::ParsedFormula<> NeumannFormula("(sin(pi*x))");

and a set of unsigned intergers identifying the edge with attribute 1, on which the Neumann boundary condition is prescribed, is created.

concepts::Set<uint> attrNeumann(static_cast<uint> (1));

Now the same is done for the Robin data. The two formulas are defined

concepts::ParsedFormula<> RobinFormula1("(-y)");
concepts::ParsedFormula<> RobinFormula2("(-sin(pi/2*y))");

and a set of unsigned intergers identifying the edge with attributes 2, 3 and 4, on which the Robin boundary condition is prescribed, is created.

concepts::Set<uint> attrRobin(static_cast<uint> (2));
attrRobin.insert(static_cast<uint> (3));
attrRobin.insert(static_cast<uint> (4));

In this example there are no homogeneous Dirichlet boundary conditions and thus, the space can be built without passing any boundary data. 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.

spc.rebuild();
graphics::drawMeshEPS(spc, "space.eps",100,1.0,1);

Now the Neumann and Robin trace spaces are built using the 2D space spc and the sets attrNeumann and attrRobin of edges with Neumann boundary condition and Robin boundary condition respectively.

hp2D::TraceSpace tspcNeumann(spc, attrNeumann);
hp2D::TraceSpace tspcRobin(spc, attrRobin);

The right hand side is computed. In our case, the vector rhs only contains the integrals of the Neumann and Robin trace as the source term is zero.

hp1D::Riesz<Real> lformNeumann(NeumannFormula);
concepts::Vector<Real> rhsNeumann(tspcNeumann, lformNeumann);
hp1D::Riesz<Real> lformRobin(RobinFormula2);
concepts::Vector<Real> rhsRobin(tspcRobin, lformRobin);
concepts::Vector<Real> rhs = rhsNeumann + rhsRobin;
std::cout << std::endl << "RHS Vector:" << std::endl << rhs << std::endl;

The system matrix is computed from the bilinear form hp2D::Laplace<Real> and the bilinear form hp2D::Identity<Real>.

A.compress();
M.compress();

Furthermore, the integrals of the Robin trace are added.

hp1D::Identity<Real> bilformRobin(RobinFormula1);
concepts::SparseMatrix<Real> M1D(tspcRobin, bilformRobin);
M1D.compress();
std::cout << std::endl << "System Matrix:" << std::endl << A << std::endl;

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

concepts::CG<Real> solver(A, 1e-6, 200);
solver(rhs, sol);
std::cout << std::endl << "Solver:" << std::endl << solver << std::endl;
std::cout << std::endl << "Solution:" << std::endl << sol << std::endl;

and plot the solution using graphics::MatlabGraphics. To this end the shape functions are computed on equidistant points using the trapezoidal quadrature rule.

spc.recomputeShapefunctions();
graphics::MatlabGraphics(spc, "Robin", sol);

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

}
std::cout << e << std::endl;
return 1;
}

# Results

Output of the program:

Mesh: Import2dMesh(ncell = 1)
RHS Vector:
Vector(169, [ 1.549354e-02, 1.678744e-01, 0.000000e+00, -9.444769e-02, 1.570060e-02, 5.980456e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -8.097468e-03, -3.200783e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 2.374103e-01, 0.000000e+00, 3.790460e-02, 2.477186e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -1.745166e-01, 0.000000e+00, 0.000000e+00, -2.305964e-02, -2.713493e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.678744e-01, 0.000000e+00, 3.790460e-02, -2.477186e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.549354e-02, -9.444769e-02, 1.570060e-02, -5.980456e-03, -8.097468e-03, -3.200783e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -1.745166e-01, 0.000000e+00, -2.305964e-02, -2.713493e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.280169e-01, -3.451119e-02, -1.813098e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.484019e-01, -2.500000e-01, -4.070872e-02, -6.366754e-04, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.500000e-01, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.280169e-01, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -3.451119e-02, -1.813098e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.500000e-01, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.484019e-01, -4.166667e-02, 1.734723e-18, -4.070872e-02, -6.366754e-04, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00])
System Matrix:
SparseMatrix(169x169, HashedSparseMatrix: 2825 (9.891110e+00%) entries bound.)
Solver:
CG(solves SparseMatrix(169x169, HashedSparseMatrix: 2825 (9.891110e+00%) entries bound.), eps = 2.020060e-13, it = 80, relres = 0)
Solution:
Vector(169, [-2.358621e-01, -1.292236e-01, -2.897387e-01, -3.566354e-01, -4.559784e-02, 1.946268e-02, -1.205115e-02, 1.041007e-03, 7.643942e-03, 1.495624e-03, 6.444425e-02, -1.673150e-02, 1.587649e-03, 1.516131e-02, -1.737078e-02, 2.614680e-04, -6.104150e-02, -2.594241e-01, 6.267788e-02, 3.415305e-03, -4.265880e-02, 3.940668e-03, 2.834395e-02, -1.218784e-03, -3.113244e-02, 3.435106e-03, 4.747729e-04, -5.287489e-04, -4.024240e-01, -4.285121e-01, -1.755820e-02, 8.663324e-04, 2.607719e-02, -8.180827e-06, -1.171876e-02, 3.220957e-04, -4.655572e-03, 1.149243e-03, -7.664004e-04, 4.331508e-05, -5.093568e-01, 3.010054e-02, -9.187840e-04, -1.329112e-02, -1.955238e-03, 1.487914e-02, -1.465611e-03, -2.471166e-03, 6.547792e-04, -1.292236e-01, -2.897387e-01, 6.267788e-02, -3.415305e-03, -1.205113e-02, 1.040996e-03, 2.834395e-02, -1.218785e-03, -3.113244e-02, 3.435105e-03, -4.747728e-04, 5.287486e-04, -2.358622e-01, -3.566353e-01, -4.559785e-02, -1.946269e-02, 6.444424e-02, -1.673150e-02, 7.643955e-03, 1.495610e-03, 1.587650e-03, 1.516131e-02, 1.737078e-02, -2.614670e-04, -5.093567e-01, -4.285121e-01, -1.329111e-02, -1.955238e-03, 3.010056e-02, -9.187924e-04, -1.171876e-02, 3.220940e-04, 1.487914e-02, -1.465611e-03, 2.471168e-03, -6.547794e-04, 2.607719e-02, -8.184208e-06, -4.655572e-03, 1.149243e-03, 7.664006e-04, -4.331492e-05, -5.422842e-01, -5.175790e-01, -1.233271e-02, 2.673123e-04, 2.503141e-02, 2.062220e-04, -1.092907e-02, -4.406529e-04, 1.897520e-03, 1.259827e-04, 3.070159e-04, -2.411602e-05, -6.223888e-01, -2.309603e-02, 9.144676e-05, 3.166487e-02, -1.023908e-03, 7.914136e-03, -6.405805e-05, 5.635018e-04, 6.261481e-05, -6.986328e-01, -6.372162e-01, -1.000691e-02, 2.816808e-03, 1.919979e-02, 1.070780e-03, -5.661673e-03, -7.812417e-04, 4.554296e-03, -2.152115e-03, 2.307448e-03, 7.109204e-06, -6.168267e-01, 2.056493e-02, -1.138131e-04, -4.696031e-03, -6.020063e-04, 8.646241e-04, -2.627329e-04, -8.207371e-05, -8.849003e-05, -5.422842e-01, -6.223888e-01, -1.233271e-02, 2.673136e-04, 3.166486e-02, -1.023899e-03, -2.309604e-02, 9.144808e-05, 7.914135e-03, -6.405777e-05, -5.635004e-04, -6.261480e-05, 2.503142e-02, 2.062226e-04, 1.897520e-03, 1.259829e-04, -3.070157e-04, 2.411587e-05, -6.372162e-01, 2.056493e-02, 1.138135e-04, -5.661684e-03, -7.812348e-04, 8.646239e-04, -2.627334e-04, 8.207376e-05, 8.849023e-05, -6.986327e-01, 1.919980e-02, -1.070775e-03, -1.000690e-02, 2.816814e-03, 4.554296e-03, -2.152114e-03, -2.307448e-03, -7.109961e-06])

Matlab plot of the solution:

@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 "hp1D.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);
// ** inhomogeneous Neumann boundary conditions at bottom using trace space **
concepts::ParsedFormula<> NeumannFormula("(sin(pi*x))");
concepts::Set<uint> attrNeumann(static_cast<uint> (1));
// ** Robin boundary conditions at right, top and left edge using trace space **
concepts::ParsedFormula<> RobinFormula1("(-y)");
concepts::ParsedFormula<> RobinFormula2("(-sin(pi/2*y))");
concepts::Set<uint> attrRobin(static_cast<uint> (2));
attrRobin.insert(static_cast<uint> (3));
attrRobin.insert(static_cast<uint> (4));
// ** space **
spc.rebuild();
graphics::drawMeshEPS(spc, "space.eps",100,1.0,1);
// ** build trace space of Neumann and Robin edges
hp2D::TraceSpace tspcNeumann(spc, attrNeumann);
hp2D::TraceSpace tspcRobin(spc, attrRobin);
// ** right hand side with added inhomogeneous Neumann BC **
hp1D::Riesz<Real> lformNeumann(NeumannFormula);
concepts::Vector<Real> rhsNeumann(tspcNeumann, lformNeumann);
hp1D::Riesz<Real> lformRobin(RobinFormula2);
concepts::Vector<Real> rhsRobin(tspcRobin, lformRobin);
concepts::Vector<Real> rhs = rhsNeumann + rhsRobin;
std::cout << std::endl << "RHS Vector:" << std::endl << rhs << std::endl;
// ** system matrix with added Robin term **
A.compress();
M.compress();
hp1D::Identity<Real> bilformRobin(RobinFormula1);
concepts::SparseMatrix<Real> M1D(tspcRobin, bilformRobin);
M1D.compress();
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;
std::cout << std::endl << "Solution:" << std::endl << sol << std::endl;
// ** plot solution **
graphics::MatlabGraphics(spc, "Robin", sol);
}
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
A function class to calculate element matrices for the mass matrix.
Definition: bilinearForm.hh:64
Builds the trace space of an FE space.
Definition: traces.hh:52
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
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 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
Linear form on edges in nD.
Definition: linearForm.hh:67
void rebuild(bool sameIndices=false)
Rebuilds the mesh and the elements due to adjustment orders.
A function class to calculate element matrices for the Laplacian.
Definition: bf_laplace.hh:91
virtual void recomputeShapefunctions()
Recompute shape functions, e.g.
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