# Introduction

One dimensional FEM 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 exact solution is u = -3 cosh(x)/cosh(1) + x^2 + 2.

The relative energy error of the FE solution is computed at the end of the program

# Commented Program

System include files for input and output operations.

#include <iostream>

Concepts include files.

#include "basics.hh"
#include "geometry.hh"
#include "space.hh"
#include "toolbox.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "linearFEM.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"

## Linear Form

The linear form is used to compute the element contributions to the load vector. The computations are done in the application operator of this class.

class RHS : public concepts::LinearForm<Real> {
private:

As we need numerical integration, the class has a member `quad_` of type concepts::Quadrature. The template parameter 4 indicates the type of quadrature rule. It contains the quadrature points and weights for given number of quadrature points (one in our case).

The constructor takes as argument the number of quadrature points and initialises the quadrature object `quad_`.

The application operator is called for every element `elm` and is responsible for computing the contributions of this element to the global load vector. The local load vector is stored in `em`.

virtual void operator()(const concepts::Element<Real>& elm,

Resize the element load vector to the correct size and set all entries to zero. This way, the computed values can be added into `em`.

const uint n = elm.T().m();
em.resize(n, 1);
em.zeros();

Check if the general element `elm` really is an element of type linearFEM::Line which can be handled by our routine. This is achieved by trying to cast `elm` to linearFEM::Line. If the cast succeeds (ie. `elm` is really of type linearFEM::Line), then `e` contains a pointer to this element. If the cast fails, `e` is equal to 0 and the assertion fails too. If the assertions fails, an error message is displayed (indicating the name of the file, the line number in the file, the name of the routine where the error happened and the check which failed) and the execution of the program stops.

const linearFEM::Line* e = dynamic_cast<const linearFEM::Line*>(&elm);

The element of type linearFEM::Line has a member `cell` which returns a reference to a concepts::Edge1d and this class is able to compute the size of the cell (on which the element is based). The size of the cell is needed to compute the Jacobian in the numerical integration.

const Real h = e->cell().size();

Next, we get pointers to the first elements of the arrays of quadrature weights and points. These data were computed in the constructor (which initialised the member `quad_`).

Now, the computational loops of this routine start.

There are two nested loops here: the outer loop enumerates the shape funtions of this element and the inner loop enumerates the quadrature points.

for (uint i = 0; i < n; ++i) { // loop over shape functions

`m` is a reference to the entry of the local load vector which corresponds to the current shape function.

Real& m = em(i,0);
for (uint k = 0; k < quad_.n(); ++k) { // loop over quadrature points

Now, all data is known to add to `m:`

• `w`[k] is the quadrature weight associated to the current quadrature point `q`[k].
m += w[k]
• `h` is the determinant of the Jacobian of the element map which maps the reference coordinates (0,1) to the physical coordinates of this element. `1/2` is the Jacobian of the map from (-1,1) to (0,1). The integration takes place on (-1,1) (as described in concepts::Quadrature).
* h / 2.0
• `e->shapefct(i, q[k])` is the value of the shape function `i` at the point `q`[k]. As you can see in linearFEM::Line, it takes values in (-1,1).
* e->shapefct(i, q[k])
• The last line `pow(e->cell().chi( (1.0+q[k])/2.0 ), 2.0)` computes the value of the right hand side funtion f = x2. `(1.0+q[k])/2.0` maps a point in (-1,1) to (0,1) and the function `pow` from the standard math library computes ab.
* pow(e->cell().chi( (1.0+q[k])/2.0 ), 2.0);
} // for k
} // for i
}
};

## Main Program

Set up the mesh. The values -1, 0 and 1 give the coordinates of the left, middle and right points in the mesh.

int main(int argc, char** argv) {
try {
#ifdef HAS_MPI
// Initialisation is necessary if MPI is installed
MPI_Init(&argc, &argv);
#endif
Line msh(-1, 0, 1);
std::cout << "Mesh: " << msh << std::endl;

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 `Boundary` holds a type (Dirichlet or Neumann) and (if necessary) a function or value. These objects can be added to `BoundaryConditions` which manages them.

std::cout << "Boundary conditions: " << bc << std::endl;

Set up the space according to the current level. The space takes the mesh, the level of refinement and the boundary conditions as input. If there are Dirichlet boundary conditions somewhere, the respective degrees of freedom are not created.

linearFEM::Linear1d spc(msh, 1, &bc);
std::cout << "Space: " << spc << std::endl;

Compute the load vector. This is done by creating an object of type `RHS` which is a linear form. This linear form is then given to the vector. The vector fills its entries by evaluating the linear form on every element of the space and assembling the local load vectors according to the assembling information contained in the elements.

RHS rhs_lf;
concepts::Vector<Real> rhs(spc, rhs_lf);

Compute the stiffness matrix and the mass matrix. As above, a bilinear form is created and given to a matrix which the assembles the mass and stiffness matrix respectively.

concepts::SparseMatrix<Real> stiff(spc, stiff_bf);
concepts::SparseMatrix<Real> mass(spc, mass_bf);

A linear combination `L` is formed from the mass and stiffness matrix. The linear system is going to be solved iteratively. This linear combination is used to initialize a CG solver together with the maximal residual and the maximal number of iterations.

concepts::LiCo<Real> L(stiff, mass, 1.0, 1.0);
concepts::CG<Real> solver(L, 1e-6, 2000);

An empty vector to hold the solution of the linear system.

The CG solver is essentially a linear operator. Therefore, one can apply a vector to it and get the solution of the linear system which was used to set up the CG solver.

try {
solver(rhs, sol);
}

The CG solver throws a exception of type `concepts::NoConvergence` which is caught here. Although the solver might not have converged, the solution might give some hints to the user. Therefore, the execution of the program is not stopped here and only a warning message is displayed.

std::cout << "solver did not converge!" << std::endl << e << std::endl;
}

Print some information about the solution process and the solution itself.

std::cout << "solver = " << solver << std::endl;
std::cout << "solution = " << sol << std::endl;

FE energy error.

const Real exEnergie = 206.0/15.0
+ 36.0*exp(2.0)/(exp(2.0)+1.0)/(exp(2.0)+1.0)*sinh(2.0)
- 72.0*exp(1.0)/(exp(2.0)+1.0)*sinh(1.0);
const Real feEnergie(sol*rhs);
std::cout << "energy error = " << fabs(exEnergie-feEnergie)/exEnergie
<< std::endl;
#ifdef HAS_MPI
MPI_Finalize();
#endif
}

Here, the exceptions are caught and some output is generated. Normally, this should not happen. We return a value different from 0 in case an exception is caught.

std::cout << e << std::endl;
return 1;
}

Sane return value.

return 0;
}

# Results

Output of the program:

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)))
Boundary conditions: BoundaryConditions(2: Boundary(DIRICHLET, (0)), 0: Boundary(FREE, (0)))
Space: Linear1d(dim = 3, nelm = 4)
solver = CG(solves LiCo(1 * SparseMatrix(3x3, HashedSparseMatrix: 7 (77.7778%) entries bound.) + 1 * SparseMatrix(3x3, HashedSparseMatrix: 7 (77.7778%) entries bound.)), eps = 9.5313e-34, it = 2, relres = 0)
solution = Vector(3, [0.0644745, 0.0642467, 0.0644745])
energy error = 0.100759
```@section complete Complete Source Code
@author Philipp Frauenfelder, 2004
```
#include <iostream>
#include "basics.hh"
#include "geometry.hh"
#include "space.hh"
#include "toolbox.hh"
#include "function.hh"
#include "operator.hh"
#include "integration.hh"
#include "linearFEM.hh"
#include "meshes.cc"
// ************************************************************* Linear Form **
class RHS : public concepts::LinearForm<Real> {
private:
public:
virtual void operator()(const concepts::Element<Real>& elm,
const uint n = elm.T().m();
em.resize(n, 1);
em.zeros();
const linearFEM::Line* e = dynamic_cast<const linearFEM::Line*>(&elm);
const Real h = e->cell().size();
for (uint i = 0; i < n; ++i) { // loop over shape functions
Real& m = em(i,0);
for (uint k = 0; k < quad_.n(); ++k) { // loop over quadrature points
m += w[k]
* h / 2.0
* e->shapefct(i, q[k])
* pow(e->cell().chi( (1.0+q[k])/2.0 ), 2.0);
} // for k
} // for i
}
};
// ************************************************************ Main Program **
int main(int argc, char** argv) {
try {
#ifdef HAS_MPI
// Initialisation is necessary if MPI is installed
MPI_Init(&argc, &argv);
#endif
Line msh(-1, 0, 1);
std::cout << "Mesh: " << msh << std::endl;
std::cout << "Boundary conditions: " << bc << std::endl;
linearFEM::Linear1d spc(msh, 1, &bc);
std::cout << "Space: " << spc << std::endl;
RHS rhs_lf;
concepts::Vector<Real> rhs(spc, rhs_lf);
concepts::SparseMatrix<Real> stiff(spc, stiff_bf);
concepts::SparseMatrix<Real> mass(spc, mass_bf);
concepts::LiCo<Real> L(stiff, mass, 1.0, 1.0);
concepts::CG<Real> solver(L, 1e-6, 2000);
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 exEnergie = 206.0/15.0
+ 36.0*exp(2.0)/(exp(2.0)+1.0)/(exp(2.0)+1.0)*sinh(2.0)
- 72.0*exp(1.0)/(exp(2.0)+1.0)*sinh(1.0);
const Real feEnergie(sol*rhs);
std::cout << "energy error = " << fabs(exEnergie-feEnergie)/exEnergie
<< std::endl;
#ifdef HAS_MPI
MPI_Finalize();
#endif
}
std::cout << e << std::endl;
return 1;
}
return 0;
}
Discrete equivalent of the Laplacian in 1D for linear FEM.
Solves a symmetric system of linear equations with conjugate gradients (CG).
Definition: cg.hh:39
Space for linear FEM in 1D.
Definition: space1D.hh:39
Exception indicating that the solver did not converge up to the desired accuracy in the given number ...
uint n() const
Returns the number of quadrature points.
Base class for exceptions.
Definition: exceptions.hh:86
void add(const Set< Attribute > &attrib, const Boundary &bcObject)
Adds a boundary condition with this attribute to the list of boundary conditions.
virtual const concepts::Edge1d & cell() const
Returns the cell of this element.
Definition: element1D.hh:53
#define conceptsAssert(cond, exc)
Assert that a certain condition is fulfilled.
Definition: exceptions.hh:394
virtual const TMatrixBase< F > & T() const =0
Returns the T matrix of the element.
Abstract class for a linear form.
Real size() const
Returns the size of the element.
const Real * weights() const
Returns a pointer into the array of the weights.
uint m() const
Returns the number of rows.
Definition: tmatrix.hh:388
Real chi(Real xi) const
The element map.
Exception class for assertions.
Definition: exceptions.hh:258
const Real * abscissas() const
Returns a pointer into the array of the abscissas.
Class to describe an element of the boundary.
Definition: boundary.hh:35
Real shapefct(const uint i, const Real xi) const
Computes the value of the ith shape function in the point xi.
virtual void operator()(const Element< G > &elm, ElementMatrix< F > &em) const =0
Computes the element contribution to the function.
virtual void resize(uint m, uint n)
Sets a new size.
Discrete equivalent for a reaction term in 1D for linear FEM.
void zeros()
Fills the matrix with zeros.
Definition: element.hh:295
Line element with linear shape function in 1D.
Definition: element1D.hh:34
Linear combination of two operators.
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