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
- Commented Program
-
Mesh
-
FEM Routine
-
Main Program
- Results
- 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>
Concepts include files.
#include "formula.hh"
#include "space.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "graphics.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:
The exact solution.
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.
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>
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.
std::cout << "Space: " << spc << std::endl;
Compute the load vector.
Compute the stiffness matrix and the mass matrix.
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.
spc.elmPairs());
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.
Solve the linear system iteratively – catch divergence errors right here.
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);
- maximal nodal error
Real maxNodalFehler = 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;
}
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.
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.
The inputParser takes care of all input parameters and the results which are gathered. At the end, they can be printed.
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.
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.
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.
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
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])
string author "(empty)"
string comment "(empty)"
string title "(empty)"
int gauss_p 1
int level 3
end
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>
#include "formula.hh"
#include "space.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "graphics.hh"
#include "meshes.cc"
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>
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;
std::cout << "Space: " << spc << std::endl;
spc.elmPairs());
spc.elmPairs());
try {
solver(rhs, sol);
}
std::cout << "solver did not converge!" << std::endl << e << std::endl;
}
std::cout << "solver = " << solver << std::endl;
const Real h = 1.0/(1<<l);
Real maxNodalFehler = 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;
}
std::stringstream filename;
filename << "solutionDG1d-" << l << ".gnuplot" << std::ends;
}
}
int main(int argc, char** argv) {
try {
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;
}
fem(input.
getInt(
"level"), input.
getInt(
"gauss_p"), output);
std::cout << inputParser << std::endl;
std::ofstream ofs("linearDG1d.gnuplot");
ofs << std::setprecision(20);
}
std::cout << e << std::endl;
return 1;
}
return 0;
}
