linearDG1d.cc

Introduction

One dimensional DGFEM in Concepts. This code features first order Ansatz functions on a (possibly) uniformly refined mesh.

The equation solved is

-u''+u = x²

on (-1,1) with homogeneous Dirichlet boundary conditions. The linear system is solved using CG. The same problem can also be solved using linear continuous FEM.

The exact solution is u = -3 cosh(x)/cosh(1) + x² + 2.

The convergence of the relative energy error and the maximal nodal error are shown together with plots of the solution.

Contents

1. Commented Program
   1. Mesh
   2. FEM Routine
   3. Main Program
2. Results
3. Complete Source Code

Commented Program

System include files. The first is for input and output operations, the second one declares auto_ptr and the third is for stringstream.

#include <iostream>
#include <memory>
#include <algorithm>
#include <sstream>
#include <unistd.h> // for command line parsing

Concepts include files.

#include "basics.hh"
#include "toolbox.hh"
#include "formula.hh"
#include "geometry.hh"
#include "space.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "graphics.hh"
#include "linearFEM.hh"
#include "linDG1D.hh"

With this using directive, we do not need to prepend concepts:: before Real everytime.

using concepts::Real;

Mesh

The mesh is taken from the meshes tutorial: Line. Here, the relevant code is included:

#include "meshes.cc"

The exact solution.

Real exact(const Real x) {
  return  -3.0*cosh(x)/cosh(1.0)+x*x+2.0;
}

FEM Routine

The only parameters of the routine fem are the maximal level of refinement and the number of quadrature points. The additional argument output is used to for the gathering the output.

void fem(const uint level, const uint gauss_p,
    concepts::InOutParameters& output) 

Set up the mesh.

{
  Line msh(-1, 0, 1);

The boundary conditions are given by assigning a boundary type to an attribute. The end points of the interval which are created above have the attribute 2.

The class concepts::Boundary holds a type (Dirichlet or Neumann) and (if necessary) a function or value. These objects can be added to concepts::BoundaryConditions which manages them.

  std::auto_ptr<concepts::BoundaryConditions>
    bc(new concepts::BoundaryConditions());
  bc->add(concepts::Attribute(2),
      concepts::Boundary(concepts::Boundary::DIRICHLET));

  std::cout << "Boundary conditions: " << *bc << std::endl;
  std::cout << "Mesh: " << msh << std::endl;

Start the loop over the different levels.

  for (uint l = 1; l < level; ++l) {
    std::cout << "l = " << l << std::endl;

And set up the space according to the current level. During this set up, the list of element pairs used later is created.

    linDG1D::Linear1d spc(msh, l);
    std::cout << "Space: " << spc << std::endl;

Compute the load vector.

    linearFEM::Riesz1d rhs_lf(concepts::ParsedFormula<>("(x^2)"), gauss_p);

Compute the stiffness matrix and the mass matrix.

    linearFEM::Laplace1d stiff_bf;
    concepts::SparseMatrix<Real> stiff(spc, stiff_bf);

    linearFEM::Identity1d mass_bf;
    concepts::SparseMatrix<Real> mass(spc, mass_bf);

Compute the boundary integrals and the stabilization integrals. These bilinear forms use a different assembly operator: it loops over the element pairs which are computed during the set up of the space. Therefore, they need to be set up differently.

    linDG1D::BoundaryInt bound_int_bf(bc.get());
    concepts::SparseMatrix<Real> bound_int(spc, bound_int_bf);
    concepts::Matrix<Real>::assembly(bound_int, bound_int_bf,
        spc.elmPairs());
    concepts::SparseMatrix<Real> bound_int_t(bound_int, true);

    linDG1D::BoundaryIntStab bound_int_stab_bf;
    concepts::SparseMatrix<Real> bound_int_stab(spc, spc);
    concepts::Matrix<Real>::assembly(bound_int_stab, bound_int_stab_bf,
        spc.elmPairs());

Finally, the different matrices are combined to form the system matrix. The parameter t can be used to chose between the SIPG (t = -1) or NIPG (t = 1) method.

    concepts::LiCo<Real> L1(stiff, mass, 1.0, 1.0);
    concepts::LiCo<Real> L2(L1, bound_int_t, 1.0, 1.0);
    const Real t = 1.0;
    concepts::LiCo<Real> L3(L2, bound_int, 1.0, t);
    concepts::LiCo<Real> L(L3, bound_int_stab, 1.0, 1.0);

Solve the linear system iteratively -- catch divergence errors right here.

    concepts::GMRes<Real> solver(L, 1e-6, 2000);
    //     concepts::SparseMatrix<Real> tmp(L);
    //     concepts::Umfpack solver(tmp);
    concepts::Vector<Real> sol(spc);
    try {
      solver(rhs, sol);
    }
    catch(concepts::NoConvergence& e) {
      std::cout << "solver did not converge!" << std::endl << e << std::endl;
    }
    std::cout << "solver = " << solver << std::endl;

Gather the results of the computation. Using output, all data can be displayed and/or saved at the end.

• mesh size h

      const Real h = 1.0/(1<<l);
      output.addArrayDouble("h", l, h);
  

• maximal nodal error

      Real x = -1.0;
      Real maxNodalFehler = 0.0;
      Real nodalFehler = 0.0;
      for (uint p = 0; p < spc.dim(); ++p) {
        nodalFehler = fabs(sol[p] - exact(x));
        if (nodalFehler > maxNodalFehler)
          maxNodalFehler = nodalFehler;
        if (p%2 == 0) x += h;
      }
      output.addArrayDouble("maxNodalFehler", l, maxNodalFehler);
  

Last but not least, the solution vector is stored in a text file suitable for Gnuplot. See the results section on how to process this file.

    std::stringstream filename;
    filename << "solutionDG1d-" << l << ".gnuplot" << std::ends;

End of the loop over the different levels.

    graphics::drawDataGnuplot(spc, filename.str(), sol);
  } // for l
}

Main Program

The main program just does some comand line parameter handling calls the main function fem. It also features the main try...catch brace to catch exceptions thrown by the program.

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

The whole program is in a big try...catch brace to catch exceptions thrown by the code.

  try {

The inputParser takes care of all input parameters and the results which are gathered. At the end, they can be printed.

    concepts::InputParser inputParser(true);
    concepts::InOutParameters& input = inputParser.inputParameters();

Default values for the variable which could be set by a command line argument: the maximal level of refinement and the number of points for the Gauss quadrature.

    input.addInt("level", 1);
    input.addInt("gauss_p", 1);

Parse the command line options (more info in input and output tutorial).

    int opt;
    while ((opt = getopt(argc, argv, "-l:p:")) != EOF)
      switch(opt) {
        case 'l': input.addInt("level", std::atoi(optarg)); break;
        case 'p': input.addInt("gauss_p", std::atoi(optarg)); break;
        default:
                  std::cout << "Usage: " << argv[0] 
                    << " [-l LEVEL] [-p GAUSS_P]" << std::endl
                    << "where" << std::endl
                    << "  LEVEL: maximal level of refinement" << std::endl
                    << "  GAUSS_P: number of points in Gauss integration"
                    << std::endl;
                  std::exit(1); break;
      }

Prepare the results part to gather the results and compute the exact FE energy and the element size.

    concepts::InOutParameters& output = inputParser.outputParameters();
    output.addArrayDouble("h");
    output.addArrayDouble("maxNodalFehler");

The values for the level and the quadrature points are fixed now, call fem.

    fem(input.getInt("level"), input.getInt("gauss_p"), output);

Print the results.

    std::cout << inputParser << std::endl;

In order to be able to plot the numbers computed in the fem routine, it is convenient to have them in table suitable for Gnuplot.

    concepts::ResultsTable table;

Every array from output which should be present in the table has to be added to table.

    table.addMap(concepts::ResultsTable::DOUBLE, "h", output);
    table.addMap(concepts::ResultsTable::DOUBLE, "maxNodalFehler", output);

Then, a file is created with high numerical output precision and the table is printed to the file. The file is closed automatically. See the results section on how to process this file.

    std::ofstream ofs("linearDG1d.gnuplot");
    ofs << std::setprecision(20);
    table.print<concepts::ResultsTable::GNUPLOT>(ofs);
  }

Here, the exceptions are caught and some output is generated. Normally, this should not happen.

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

Sane return value

  return 0;
}

Results

Output of the program:

$ linearDG1d -l 3
Boundary conditions: BoundaryConditions(2: Boundary(DIRICHLET, (0)), 0: Boundary(FREE, (0)))
Mesh: Line(Edge1d(idx = (0, 0), cntr = Edge(Key(0), (Vertex(Key(0)), Vertex(Key(1))), Attribute(0)), vtx = [-1, 0], map =
MapEdge1d(-1, 0)), Edge1d(idx = (0, 0), cntr = Edge(Key(1), (Vertex(Key(1)), Vertex(Key(2))), Attribute(0)), vtx = [0, 1], map = MapEdge1d(0, 1)))
l = 1
Space: Linear1d(dim = 8, nelm = 4)
solver = GMRes(solves LiCo(1 * LiCo(1 * LiCo(1 * LiCo(1 * SparseMatrix(8x8, HashedSparseMatrix: 16 (25%) entries bound.) + 1 * SparseMatrix(8x8, HashedSparseMatrix: 16 (25%) entries bound.)) + 1 * SparseMatrix(8x8, HashedSparseMatrix: 28 (43.75%) entries bound.)) + -1 * SparseMatrix(8x8, HashedSparseMatrix: 40 (62.5%) entries bound.)) + 1 * SparseMatrix(8x8, HashedSparseMatrix: 40 (62.5%) entries bound.)), eps = 1.01899e-16, it = 4)
solution = Vector(8, [0.0394152, 0.11776, 0.0843423, 0.0663982, 0.0663982, 0.0843423, 0.11776, 0.0394152])
l = 2
Space: Linear1d(dim = 16, nelm = 8)
solver = GMRes(solves LiCo(1 * LiCo(1 * LiCo(1 * LiCo(1 * SparseMatrix(16x16, HashedSparseMatrix: 32 (12.5%) entries bound.) + 1 * SparseMatrix(16x16, HashedSparseMatrix: 32 (12.5%) entries bound.)) + 1 * SparseMatrix(16x16, HashedSparseMatrix: 60 (23.4375%) entries bound.)) + -1 * SparseMatrix(16x16, HashedSparseMatrix: 88 (34.375%) entries bound.)) + 1 * SparseMatrix(16x16, HashedSparseMatrix: 88 (34.375%) entries bound.)), eps = 2.35627e-16, it = 8)
solution = Vector(16, [0.0135647, 0.0676432, 0.0592614, 0.0666052, 0.0655731, 0.0639907, 0.0628846, 0.0608961, 0.0608961,
0.0628846, 0.0639907, 0.0655731, 0.0666052, 0.0592614, 0.0676432, 0.0135647])
/* time, date:  Mon Aug  4 15:33:40 2003
 username:    pfrauenf
 sysname:     Linux
 nodename:    herodot
 release:     2.4.20
 version:     2.4.20
 machine:     i686
 */
string  author  "(empty)"
string  comment "(empty)"
string  title   "(empty)"
int     gauss_p 1
int     level   3
end
// output starts here
array double h {
        1       0.5
        2       0.25
}
array double maxNodalFehler {
        1       0.060049
        2       0.0222183
} 

Error diagram. The following diagram shows the convergence rates of the method symmetric (t = -1) and the unsymmetric (t = 1) method.

The diagram above was created using the following commands:

set logscale x
set logscale y
set grid
set grid mxtics
set grid mytics
set key top left
set xlabel 'h'
set ylabel 'error'
set data style linespoints
plot 'linearDG1d_t1.gnuplot' using 2:3 title 'maximal nodal error, t=1', \
  'linearDG1d_t-1.gnuplot' using 2:3 title 'maximal nodal error, t=-1' 

Complete Source Code

Author:

Philipp Frauenfelder, 2004

#include <iostream>
#include <memory>
#include <algorithm>
#include <sstream>
#include <unistd.h> // for command line parsing

#include "basics.hh"
#include "toolbox.hh"
#include "formula.hh"
#include "geometry.hh"
#include "space.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "graphics.hh"
#include "linearFEM.hh"
#include "linDG1D.hh"

using concepts::Real;

#include "meshes.cc"

Real exact(const Real x) {
  return  -3.0*cosh(x)/cosh(1.0)+x*x+2.0;
}

void fem(const uint level, const uint gauss_p,
    concepts::InOutParameters& output) 
{
  Line msh(-1, 0, 1);

  std::auto_ptr<concepts::BoundaryConditions>
    bc(new concepts::BoundaryConditions());
  bc->add(concepts::Attribute(2),
      concepts::Boundary(concepts::Boundary::DIRICHLET));

  std::cout << "Boundary conditions: " << *bc << std::endl;
  std::cout << "Mesh: " << msh << std::endl;

  for (uint l = 1; l < level; ++l) {
    std::cout << "l = " << l << std::endl;

    linDG1D::Linear1d spc(msh, l);
    std::cout << "Space: " << spc << std::endl;

    linearFEM::Riesz1d rhs_lf(concepts::ParsedFormula<>("(x^2)"), gauss_p);
    concepts::Vector<Real> rhs(spc, rhs_lf);

    linearFEM::Laplace1d stiff_bf;
    concepts::SparseMatrix<Real> stiff(spc, stiff_bf);

    linearFEM::Identity1d mass_bf;
    concepts::SparseMatrix<Real> mass(spc, mass_bf);

    linDG1D::BoundaryInt bound_int_bf(bc.get());
    concepts::SparseMatrix<Real> bound_int(spc, bound_int_bf);
    concepts::Matrix<Real>::assembly(bound_int, bound_int_bf,
        spc.elmPairs());
    concepts::SparseMatrix<Real> bound_int_t(bound_int, true);

    linDG1D::BoundaryIntStab bound_int_stab_bf;
    concepts::SparseMatrix<Real> bound_int_stab(spc, spc);
    concepts::Matrix<Real>::assembly(bound_int_stab, bound_int_stab_bf,
        spc.elmPairs());

    concepts::LiCo<Real> L1(stiff, mass, 1.0, 1.0);
    concepts::LiCo<Real> L2(L1, bound_int_t, 1.0, 1.0);
    const Real t = 1.0;
    concepts::LiCo<Real> L3(L2, bound_int, 1.0, t);
    concepts::LiCo<Real> L(L3, bound_int_stab, 1.0, 1.0);

    concepts::GMRes<Real> solver(L, 1e-6, 2000);
    //     concepts::SparseMatrix<Real> tmp(L);
    //     concepts::Umfpack solver(tmp);
    concepts::Vector<Real> sol(spc);
    try {
      solver(rhs, sol);
    }
    catch(concepts::NoConvergence& e) {
      std::cout << "solver did not converge!" << std::endl << e << std::endl;
    }
    std::cout << "solver = " << solver << std::endl;
    //     std::cout << "solution = " << sol << std::endl;

    const Real h = 1.0/(1<<l);
    output.addArrayDouble("h", l, h);

    Real x = -1.0;
    Real maxNodalFehler = 0.0;
    Real nodalFehler = 0.0;
    for (uint p = 0; p < spc.dim(); ++p) {
      nodalFehler = fabs(sol[p] - exact(x));
      if (nodalFehler > maxNodalFehler)
        maxNodalFehler = nodalFehler;
      if (p%2 == 0) x += h;
    }
    output.addArrayDouble("maxNodalFehler", l, maxNodalFehler);

    std::stringstream filename;
    filename << "solutionDG1d-" << l << ".gnuplot" << std::ends;
    graphics::drawDataGnuplot(spc, filename.str(), sol);
  } // for l
}

// ************************************************************ Main Program **

int main(int argc, char** argv) {
  try {
    concepts::InputParser inputParser(true);
    concepts::InOutParameters& input = inputParser.inputParameters();

    input.addInt("level", 1);
    input.addInt("gauss_p", 1);

    int opt;
    while ((opt = getopt(argc, argv, "-l:p:")) != EOF)
      switch(opt) {
        case 'l': input.addInt("level", std::atoi(optarg)); break;
        case 'p': input.addInt("gauss_p", std::atoi(optarg)); break;
        default:
                  std::cout << "Usage: " << argv[0] 
                    << " [-l LEVEL] [-p GAUSS_P]" << std::endl
                    << "where" << std::endl
                    << "  LEVEL: maximal level of refinement" << std::endl
                    << "  GAUSS_P: number of points in Gauss integration"
                    << std::endl;
                  std::exit(1); break;
      }

    concepts::InOutParameters& output = inputParser.outputParameters();
    output.addArrayDouble("h");
    output.addArrayDouble("maxNodalFehler");

    fem(input.getInt("level"), input.getInt("gauss_p"), output);

    std::cout << inputParser << std::endl;

    concepts::ResultsTable table;
    table.addMap(concepts::ResultsTable::DOUBLE, "h", output);
    table.addMap(concepts::ResultsTable::DOUBLE, "maxNodalFehler", output);
    std::ofstream ofs("linearDG1d.gnuplot");
    ofs << std::setprecision(20);
    table.print<concepts::ResultsTable::GNUPLOT>(ofs);
  }
  catch(concepts::ExceptionBase& e) {
    std::cout << e << std::endl;
    return 1;
  }

  return 0;
}