Numerical C++ Library Concepts

TWiki Guest (guest)

Contents

Introduction
Variational Formulation
Commented Program
Results
Complete Source Code

Introduction

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

The equation solved is the reaction-diffusion equation

$- \Delta u + u = 0$

in the unit square $\Omega = (0,1)^2$ with inhomogeneous Neumann boundary condition

$\nabla u\cdot \underline{n} = g(x,y), \qquad (x,y)\in\Gamma_{\rm N}=\{(x,y)\in\partial\Omega\mid y=0\},$

and homogeneous Dirichlet boundary conditions

$u = 0, \qquad (x,y)\in\Gamma_{\rm D}=\partial\Omega\setminus\Gamma_{\rm N},$

where $g(x,y) = \sin \pi x$ .

Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find $u\in X=\{v\in H^1(\Omega)\mid \left.v\right|_{\Gamma_{\rm D}}=0\}$ such that

$\int_{\Omega} \nabla u \cdot \nabla v \;{\rm d}x + \int_{\Omega} u v \;{\rm d}x = \int_{\Gamma_{\rm N}}gv\;{\rm d}x \qquad\forall v\in X.$

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

using concepts::Real;

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);

The homogeneous Dirichlet boundary conditions are set.

    concepts::BoundaryConditions bc;
    for (uint i = 2; i <= 4; ++i)
      bc.add(concepts::Attribute(i), concepts::Boundary(concepts::Boundary::DIRICHLET));

The formula of the inhomogeneous Neumann data is defined

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

and a set of unsigned intergers identifying the edges with Neumann boundary conditions is created.

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

Note that in our example Neumann boundary conditions are only prescribed at the bottom boundary with attribute 1.

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.

    hp2D::hpAdaptiveSpaceH1 spc(msh, 2, 3, &bc);
    spc.rebuild();
    graphics::drawMeshEPS(spc, "space.eps",100,1.0,1);

Now the Neumann trace space is built using the 2D space spc and the set attrNeumann of edges with Neumann boundary condition.

hp2D::TraceSpace tspcNeumann(spc, attrNeumann);

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

    hp1D::Riesz<Real> lformNeumann(NeumannFormula);
    concepts::Vector<Real> rhs(tspcNeumann, lformNeumann);
    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>.

    hp2D::Laplace<Real> la;
    concepts::SparseMatrix<Real> A(spc, la);
    A.compress();
    hp2D::Identity<Real> id;
    concepts::SparseMatrix<Real> M(spc, id);
    M.compress();
    M.addInto(A, 1.0);
    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);
    concepts::Vector<Real> sol(spc);
    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.

    hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, 8);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "inhomNeumann", sol);

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

  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }

Results

Output of the program:

Mesh: Import2dMesh(ncell = 1)
 
RHS Vector:
Vector(132, [ 1.678744e-01,  0.000000e+00,  1.570060e-02,  5.980456e-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,  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,  0.000000e+00,  0.000000e+00,  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.570060e-02, -5.980456e-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,  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,  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])
 
System Matrix:
SparseMatrix(132x132, HashedSparseMatrix: 1794 (1.029614e+01%) entries bound.)
 
Solver:
CG(solves SparseMatrix(132x132, HashedSparseMatrix: 1794 (1.029614e+01%) entries bound.), eps = 7.757995e-13, it = 36, relres = 0)
 
Solution:
Vector(132, [ 2.138848e-01,  9.326628e-02,  3.516844e-02,  5.584007e-03, -4.883629e-02,  3.356455e-03,  1.532876e-02,  2.433166e-03, -8.064913e-03,  5.630397e-04, -1.288910e-03,  9.407855e-05,  3.024788e-01,  1.318984e-01,  8.490413e-02,  2.312972e-03, -6.906494e-02,  4.746742e-03,  3.700689e-02, -1.007849e-03, -1.947042e-02,  1.359298e-03, -5.338839e-04,  3.896852e-05,  5.610815e-02,  3.967446e-02, -2.990543e-02,  2.106316e-03,  1.573143e-02,  4.281191e-04, -2.114634e-02, -1.489391e-03, -8.355282e-03,  6.177991e-04, -2.262989e-04,  1.911231e-05,  6.516172e-03,  1.033572e-03, -3.460871e-03,  2.559008e-04, -5.463338e-04,  4.614120e-05,  2.138848e-01,  9.326629e-02,  8.490413e-02, -2.312971e-03, -4.883628e-02,  3.356452e-03,  3.700689e-02, -1.007851e-03, -1.947042e-02,  1.359298e-03,  5.338840e-04, -3.896865e-05,  3.516844e-02, -5.584008e-03,  1.532875e-02,  2.433164e-03, -8.064912e-03,  5.630394e-04,  1.288910e-03, -9.407842e-05,  3.967446e-02,  6.516172e-03,  1.033571e-03, -2.114633e-02, -1.489390e-03, -3.460871e-03,  2.559007e-04,  5.463338e-04, -4.614119e-05,  1.573143e-02,  4.281196e-04, -8.355282e-03,  6.177991e-04,  2.262988e-04, -1.911231e-05,  1.459316e-02,  2.063785e-02, -8.635381e-03,  6.999117e-04,  5.786236e-03, -1.585277e-04, -1.221227e-02, -9.898247e-04, -3.384507e-03,  2.667914e-04,  9.347732e-05, -9.106996e-06,  2.396737e-03,  3.827195e-04, -1.401909e-03,  1.105086e-04,  2.256743e-04, -2.198624e-05, -2.321856e-03, -4.072598e-04, -3.763560e-04,  6.414338e-05,  7.463472e-05, -1.439713e-05, -3.283600e-03, -5.759523e-04, -9.086037e-04,  1.548558e-04,  3.091471e-05, -5.963484e-06,  1.459316e-02, -8.635382e-03,  6.999119e-04,  2.396737e-03,  3.827196e-04, -1.401909e-03,  1.105086e-04, -2.256742e-04,  2.198624e-05,  5.786236e-03, -1.585275e-04, -3.384507e-03,  2.667914e-04, -9.347736e-05,  9.107000e-06, -2.321856e-03, -4.072598e-04, -9.086037e-04,  1.548558e-04, -3.091471e-05,  5.963487e-06, -3.763560e-04,  6.414338e-05, -7.463471e-05,  1.439713e-05])

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 "hp1D.hh"
#include "hp2D.hh"
#include "operator.hh"
#include "space.hh"
#include "toolbox.hh"
 
using concepts::Real;
 
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);
 
    // ** homogeneous Dirichlet boundary conditions at the right, top and left edge **
    concepts::BoundaryConditions bc;
    for (uint i = 2; i <= 4; ++i)
      bc.add(concepts::Attribute(i), concepts::Boundary(concepts::Boundary::DIRICHLET));
    
    // ** inhomogeneous Neumann boundary conditions at the bottom using trace space **
    concepts::ParsedFormula<> NeumannFormula("(sin(pi*x))");
    concepts::Set<uint> attrNeumann(static_cast<uint> (1));
    
    // ** space **
    hp2D::hpAdaptiveSpaceH1 spc(msh, 2, 3, &bc);
    spc.rebuild();
    graphics::drawMeshEPS(spc, "space.eps",100,1.0,1);
    
    // ** build trace space of Neumann edge
    hp2D::TraceSpace tspcNeumann(spc, attrNeumann);
        
    // ** right hand side with added inhomogeneous Neumann BC **
    hp1D::Riesz<Real> lformNeumann(NeumannFormula);
    concepts::Vector<Real> rhs(tspcNeumann, lformNeumann);
    std::cout << std::endl << "RHS Vector:" << std::endl << rhs << std::endl;
    
    // ** left hand side matrix **    
    hp2D::Laplace<Real> la;
    concepts::SparseMatrix<Real> A(spc, la);
    A.compress();
    hp2D::Identity<Real> id;
    concepts::SparseMatrix<Real> M(spc, id);
    M.compress();
    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);
    concepts::Vector<Real> sol(spc);
    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 **  
    hp2D::IntegrableQuad::factory().get()->setTensor(concepts::TRAPEZE, true, 8);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "inhomNeumann", sol);
  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }
  return 0;
}

graphics::MatlabGraphics

Draws a picture of data in Matlab format and stores the result in a single file.

Definition: matlab.hh:114

hp2D::hpAdaptiveSpaceH1

Definition: hpAdaptiveSpaceH1.hh:63

concepts::BoundaryConditions

Boundary conditions.

Definition: boundaryConditions.hh:38

hp2D::TraceSpace

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

hp2D::IntegrableQuad::factory

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

concepts::TRAPEZE

@ TRAPEZE

Definition: defines.hh:13

concepts::ParsedFormula<>

concepts::Import2dMesh

Definition: meshImport2D.hh:103

concepts::ExceptionBase

Base class for exceptions.

Definition: exceptions.hh:86

concepts::Boundary::DIRICHLET

@ DIRICHLET

Definition: boundary.hh:38

graphics::drawMeshEPS

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.

concepts::BoundaryConditions::add

void add(const Set< Attribute > &attrib, const Boundary &bcObject)

Adds a boundary condition with this attribute to the list of boundary conditions.

concepts::Set< uint >

concepts::SparseMatrix::compress

void compress(Real threshold=EPS)

Compresses the matrix by dropping small entries.

Definition: sparseMatrix.hh:444

concepts::SparseMatrix::addInto

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.

Definition: sparseMatrix.hh:538

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

concepts::Boundary

Class to describe an element of the boundary.

Definition: boundary.hh:35

hp2D::hpAdaptiveSpace::rebuild

void rebuild(bool sameIndices=false)

Rebuilds the mesh and the elements due to adjustment orders.

concepts::Vector< Real >

A function class to calculate element matrices for the Laplacian.

Definition: bf_laplace.hh:91

hp2D::hpAdaptiveSpaceH1::recomputeShapefunctions

virtual void recomputeShapefunctions()

Recompute shape functions, e.g.

concepts::Attribute

Attributes for elements of the topology.

Definition: connector.hh:22

double Real

Type normally used for a floating point number.

Definition: typedefs.hh:36

concepts::SparseMatrix< Real >

Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich