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 = x2

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) + x2 + 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 unique_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.

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,

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::unique_ptr<concepts::BoundaryConditions>
bc->add(concepts::Attribute(2),
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.

concepts::SparseMatrix<Real> stiff(spc, stiff_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);
try {
solver(rhs, sol);
}
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.

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

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

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

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.
<img src="linearDG1d-error.png" width="640" height="480" alt="error diagram">
The diagram above was created using the following commands: @code

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"
#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,
{
Line msh(-1, 0, 1);
std::unique_ptr<concepts::BoundaryConditions>
bc->add(concepts::Attribute(2),
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);
concepts::SparseMatrix<Real> stiff(spc, stiff_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);
try {
solver(rhs, sol);
}
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;
table.addMap(concepts::ResultsTable::DOUBLE, "maxNodalFehler", output);
std::ofstream ofs("linearDG1d.gnuplot");
ofs << std::setprecision(20);
}
std::cout << e << std::endl;
return 1;
}
return 0;
}
Discrete equivalent of the Laplacian in 1D for linear FEM.
Exception indicating that the solver did not converge up to the desired accuracy in the given number ...
Base class for exceptions.
Definition: exceptions.hh:86
Holds parameters in hashes.
Definition: inputOutput.hh:75
Boundary integral for the DG FEM in 1D.
Definition: bilinearForm.hh:27
void addArrayDouble(const char *array, const bool newArray=false)
Creates an empty array for doubles if necessary.
Linear form form 1D linear FEM.
Definition: linearForm1D.hh:26
void drawDataGnuplot(T &msh_spc, std::string filename, const S &sol, const concepts::ElementFunction< Real > *fun=0)
Creates a data file for viewing the data with Gnuplot using DataGnuplotCell.
Definition: dataGnuplot.hh:109
int getInt(const char *name, const int value=INT_MAX) const
Returns an int from the hash of ints.
Class to describe an element of the boundary.
Definition: boundary.hh:35
Parses an input file and extracts parameters.
Definition: inputOutput.hh:443
void addMap(enum mapTypes type, const char *name, const InOutParameters &holder)
void addInt(const char *name, const int value)
Adds an int to the hash of ints.
Space for linear DG FEM in 1D.
Definition: space.hh:34
Stabilizing boundary integral for the DG FEM in 1D.
Definition: bilinearForm.hh:49
Discrete equivalent for a reaction term in 1D for linear FEM.
Organizes the results in the hashes from InOutParameters in a nice table.
Definition: resultsTable.hh:23
Solves a system of linear equations with general minimal residuals (GMRes).
Definition: gmres.hh:24
static void assembly(Matrix< F > &dest, const Space< G > &spc, const BilinearForm< F, G > &bf, const Real threshold=0.0)
Assembly operator for dest using the bilinear form bf.
Linear combination of two operators.
void print(std::ostream &os) const
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: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich