Home | Doxygen Documentation | Tutorials | Developer Tools (restricted)
inhomDirichletBCs.cc
  1. Introduction
  2. Variational Formulation
  3. Commented Program
  4. Results
  5. Complete Source Code

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

\[ - \Delta u + u = f \]

in the unit square $ \Omega = (0,1)^2 $ with source term

\[ f(x,y) = -\exp\left(-(x-1/2)^2-(y-1/2)^2\right) \]

and Dirichlet boundary conditions

\[ u = \begin{cases} \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}. \end{cases} \]

A suitable Dirichlet lift is

\[ g = \sin \pi x \; \cos \frac{\pi}{2} y. \]

By setting $ u = \tilde{u} + g $ the reaction-diffusion equation simplifies to a problem in $ \tilde{u} $

\[ - \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. \]

Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find $ \tilde{u}\in H^1_0(\Omega) $ 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). \]

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 "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 Dirichlet lift and its gradient are defined. 

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

The problem for $ \tilde{u} $ has homogeneous Dirichlet boundary conditions at all four edges. 

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

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

The right hand side is computed, containing the integrals of the source term, the Dirichlet lift and its gradient. 

    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;

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;

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 spc, into an element formula, that can be evaluated at each point in every cell. 

    concepts::ElementFormulaVector<1> solFormula(spc, sol, hp2D::Value<Real>());
    std::cout << std::endl << "Solution:" << std::endl << solFormula + DirichletLift << std::endl;

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

    hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, 8);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "inhomDirichlet", solFormula + DirichletLift);
    graphics::MatlabGraphics(spc, "inhomDirichlet0", solFormula);
    graphics::MatlabGraphics(spc, "inhomDirichletLift", DirichletLift);

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(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:

inhomDirichlet.png

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"

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

    // ** 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 **
    concepts::BoundaryConditions bc;
    for (uint i = 1; i <= 4; ++i)
      bc.add(concepts::Attribute(i), concepts::Boundary(concepts::Boundary::DIRICHLET));
    
    // ** space **
    hp2D::hpAdaptiveSpaceH1 spc(msh, 2, 3, &bc);
    spc.rebuild();
    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 **
    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;
    
    // ** transform solution to formula (in order to add Dirichlet lift) **
    concepts::ElementFormulaVector<1> solFormula(spc, sol, hp2D::Value<Real>());
    std::cout << std::endl << "Solution:" << std::endl << solFormula + DirichletLift << std::endl;
        
    // ** plot solution **  
    hp2D::IntegrableQuad::rule().set(concepts::TRAPEZE, true, 8);
    spc.recomputeShapefunctions();
    graphics::MatlabGraphics(spc, "inhomDirichlet", solFormula + DirichletLift);
    graphics::MatlabGraphics(spc, "inhomDirichlet0", solFormula);
    graphics::MatlabGraphics(spc, "inhomDirichletLift", DirichletLift);
  }
  catch (concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }
  return 0;
}
Home | Doxygen Documentation | Tutorials | Developer Tools (restricted)