# Tutorial 00: P1 Finite Elements for the Poisson Equation

## Problem Formulation

In this exercise, we will solve our first PDE, the *Poisson equation*:

\begin{align}\label{eq:1a}
\begin{array}{rcll}
-\Delta u & = f && \text{in}\ \Omega, \\
u & = g && \text{on}\ \partial\Omega, 
\end{array}
\end{align}


where $\Omega \in \mathbb{R}^d$ is a given, polyhedral domain
(elements with curved boundaries are possible in DUNE but will not be considered
here). This problem is one of the basic equations of mathematical physics
which describes gravitational and electric potential as well as stationary heat or
groundwater flow. Poisson's equation is an instance of an *elliptic partial differential equation*.

A function $u \in C^2(\Omega) \cap C^0(\bar{\Omega})$ (the space of twice continuously differentiable functions in $\Omega$ which are continuous up to the boundary) is called a *strong solution* if it satisfies the above equations pointwise.

Formally, one often reduces the equations to a problem with homogeneous Dirichlet boundary conditions g ≡ 0 as the starting point for theoretical considerations. We will deliberately not do this here, as it is also not done in the computer implementation of the method.

### Weak Formulation
<a id='weakFormulation'> </a>

Existence, uniqueness and stability of solutions, i.e. well-posedness in the sense of
Hadamard, is easier to prove for so-called weak solutions. As the weak formulation
is also the basis of the finite element method we explain it here.

As a start, suppose $u$ is a strong solution of the equations above and take 
*any* function
$v\in C^1(\Omega)\cap C^0(\bar\Omega)$ with $v=0$ on $\partial\Omega$, then we
have by integration by parts:
\begin{equation}
\int_\Omega (-\Delta u) v \,dx = 
\int_\Omega \nabla u \cdot \nabla v \,dx= \int_\Omega fv \,dx.
\end{equation}
Observe, that the boundary integral $\int_{\partial\Omega} \nabla u\cdot n v\,dx$ vanishes
due to the fact that $v=0$ on $\partial\Omega$. Loosely speaking $v=0$ is a consequence
of the Dirichlet boundary condition $u=g$ on $\partial\Omega$.

Introducing the abbreviations
\begin{equation}
a(u,v) = \int_\Omega \nabla u \cdot \nabla v \,dx, \qquad l(v) = \int_\Omega fv \,dx
\end{equation}
one can on the other hand ask the question: Is there a class,
more specific a vector space, of functions $V$ with $V_g=\{v\in V : 
\text{$v=g$ on $\partial\Omega$}\}$ and $V_0=\{v\in V : 
\text{$v=0$ on $\partial\Omega$}\}$ such that
the problem
\begin{equation} 
u \in V_g :\qquad a(u,v) = l(v) \qquad \forall v\in V_0 \label{eq:weakform}
\end{equation}
has a unique solution.

The answer is yes and in particular one can prove the following:

1. With $V=H^1(\Omega)$, the Sobolev space of functions with square integrable
weak derivatives, the problem \eqref{eq:weakform} has a unique solution provided
the bilinear form $a : V\times V \to \mathbb{R}$ is continuous and coercive on the 
subspace $V_0\subset V$ and the linear form $l: V \to \mathbb{R}$ is also continuous. 
Coercivity on $V_0$ follows from Friedrich's inequality and for the continuity of the right hand
side functional $f\in L^2(\Omega)$ is a sufficient condition.
2. If in addition $u\in C^2(\Omega)\cap C^0(\bar\Omega)$, then the solutions
of \eqref{eq:weakform} and \eqref{eq:1a} coincide.

We call \eqref{eq:weakform} the *weak formulation* of \eqref{eq:1a}.
As it has a unique solution under more general conditions than \eqref{eq:1a},
e.g. for discontinuous right hand side functions $f$, it can be considered a generalization
of the problem.

## The Finite Element Method

The (conforming) finite element method, in a nutshell, is based on the weak formulation where
the function space $V$ is replaced by a subspace $V_h\subset V$ 
which is *finite dimensional*. Here, the subscript $h$ relates to the dimension of the
function space. One major part of the finite element method is to construct such so-called *finite element spaces*.
Typically, finite element spaces consist of piecewise polynomial functions.
We consider one particular choice, the space of linear functions on simplicial elements.
For a more detailed description of the underlying finite element method, please refer to `tutorial00.pdf`. Additionally, there are many text books about the finite element method, which are listed in the [further reading section](#FurtherReading) at the end of this tutorial.

## Implementation

The code presented in this tutorial solves the elliptic equation \eqref{eq:1a} with Dirichlet boundary conditions, where
\begin{align}
    f(x) = -2d \quad\text{and}\quad  g(x) = \sum_{i=1}^d (x)_i^2.
\end{align}

*Note that throughout this section there are 'show solution buttons' included. 
They only become relevant for the exercises in [section 2](#Exercises), as these exercises require some modification within the code.*

First, we include all the Dune sources through one convenience header:


```c++
#include<dune/jupyter.hh>
```

The simulation code written below will depend on a number of runtime parameters that are provided in a configuration file `tutorial00.ini`. You find that file in the `notebooks` directory. Dune uses a nested data structured called `ParameterTree` for these configurations. In the following we parse the configuration, output it to the notebook and alter it from within the notebook:


```c++
Dune::ParameterTree ptree;
Dune::ParameterTreeParser::readINITree("tutorial00.ini", ptree);
```


```c++
ptree
```


```c++
ptree["grid.refinement"] = "0";
```


```c++
ptree
```

We are using a two-dimensional, unstructured mesh for this simulation. The mesh resolves the unit square $\Omega = [0,1]^2$. The file is specified in the `unitsquare.msh` file that was generated using Gmsh.


```c++
const int dim = 2;
using Grid = Dune::UGGrid<dim>;
std::string filename = ptree.get("grid.twod.filename", "unitsquare.msh");
Dune::GridFactory<Grid> factory;
Dune::GmshReader<Grid>::read(factory,filename,true,true);
std::unique_ptr<Grid> grid(factory.createGrid());
grid->globalRefine(ptree.get<int>("grid.refinement"));
auto gv = grid->leafGridView();
using GV = decltype(gv);
```


```c++
grid
```

With Dune making heavy use of C++ templates, the underlying floating point type of the simulation can be chosen quite freely. Looking for a solution $u: \Omega\rightarrow\mathbb{R}$, the *domain field* type is the floating point type used to realize $\Omega$ where as the *range field* type is used to realize $\mathbb{R}$. `DF` is defined by the grid implementation (usually to `double`), but we can freely choose `RF` - typically also to `double`:


```c++
using DF = Grid::ctype;
using RF = double;
```

In the following, first a *local finite element map* of type `PkLocalFiniteElementMap` is set up. The finite element map associates finite element basis functions, defined on the corresponding reference element, with each element of the mesh. Additionally, information is provided on how the global finite element space $V_h$ is constructed from the local, element-wise pieces.

The type `CON` provides a so-called constraints class, which provides a way to assemble constraints on a function space.

The type `VBE` provides a vector backend, this allows for the direct filling of different types of data into these vectors, without a copy step. This is usefull, as PDELab can use different iterative solver libraries, which provide their own data types for vectors and (sparse) matrices. In this code, DUNE's own iterative solver library ISTL provides the vector backend.

From the above described parameters, the type `GFS` from the class template `GridFunctionSpace`is defined. This combines all the given information to construct the global finite element space $V_h$ on a given grid view. In the last two lines an object of this class is instantiated and the space is given a name.


```c++
using FEM = Dune::PDELab::PkLocalFiniteElementMap<GV,DF,RF,1>;
FEM fem(gv);
using CON = Dune::PDELab::ConformingDirichletConstraints;
using VBE = Dune::PDELab::ISTL::VectorBackend<>;
using GFS = Dune::PDELab::GridFunctionSpace<GV,FEM,CON,VBE>;
GFS gfs(gv,fem);
gfs.name("P1");
```

In the next step, a variable that will later on contain the solution vector *z* is declared:


```c++
using Z = Dune::PDELab::Backend::Vector<GFS,RF>;
Z z(gfs);
```

The function g extends the Dirichlet boundary values into the inerior. This function can be used to provide e.g. the exact solution of the problem for testing purposes or an initial guess for the iterative solvers. *Note that as stated above, one exercise at the end of this notebook requires a change of this function. The solution for this part can be accessed through the "show solution" button below.*
<a id='newfandg'> </a>


```c++
auto g = Dune::PDELab::makeGridFunctionFromCallable(
    gv,
    [](const auto& x){
        RF s=0.0;
        for (std::size_t i=0; i<x.size(); i++)
            s+=x[i]*x[i];
        return s;
    }
);
Dune::PDELab::interpolate(g,gfs,z);
```

*This is the solution for exercise 1.2. Note that you need to comment out the above lines before executing the notebook with the new function g.*


```c++
//solution to exercise 1.2
/*
auto g = Dune::PDELab::makeGridFunctionFromCallable(
    gv,
    [](const auto& x){
        RF s=0.0;
        for (std::size_t i=0; i<x.size(); i++)
            s+=x[i]*x[i]*x[i];
        return s;
    }
);
Dune::PDELab::interpolate(g,gfs,z);
*/
```

In the following, the constraints for a specific grid function space are determined from the boundary condition type function `b`. They are saved into the constraints container `cc` of type `CC`. The separation of the function space $V_h$ and the constraints set $\mathcal{I}_h^{\partial \Omega}$ allows one to reuse the function space together with different constraints, e.g. for a system of partial differential equations.


```c++
using CC = typename GFS::template ConstraintsContainer<RF>::Type;
CC cc;
```


```c++
auto b = Dune::PDELab::makeBoundaryConditionFromCallable(
    gv,
    [](const auto& x){ return true; }
);
Dune::PDELab::constraints(b, gfs, cc);
std::cout << "constrained dofs=" << cc.size() << " of "
          << gfs.globalSize() << std::endl;
```

The right hand side $f$ of the partial differential equation is denoted by the function `f`. *Note that the exercise also requires a change of the function f, thus there is another solution.*


```c++
auto f = Dune::PDELab::makeGridFunctionFromCallable(
    gv,
    [](const auto& x){
        return Dune::FieldVector<RF,1>(-2.0*x.size());
    }
);;
```

*This is the solution to exercise 1.2. Again make sure only one of the cells containg the function f is executed.*


```c++
//solution to exercise 1.2
/*
auto f = Dune::PDELab::makeGridFunctionFromCallable(
    gv,
    [](const auto& x){
        RF s=0.0;
        RF t=0.0;
        for (std::size_t i=0; i<x.size(); i++){
          s+=x[i];
          t+=x[i]*x[i]*x[i];
        }
        return -6*s;
    }
);;
*/
```

*Exercise 2 requires additional changes, the solution can be found below.*
<a id='solex2'> </a>


```c++
/*
auto k = Dune::PDELab::makeGridFunctionFromCallable(
    gv, 
    [](const auto& x){
      return x[0];
    }
);;
*/
```


```c++
/*
auto a = Dune::PDELab::makeGridFunctionFromCallable(
   gv, 
   [](const auto& x) {
       RF s=0.0;
       for (std::size_t i=0; i<x.size(); i++)
         s+=x[i];
       return s;
   }
);;
*/
```

Next, we will construct the local operator. We include the local operator implementation for P1-conforming finite elements from the file `poissonp1.hh`. You find this file in the `notebooks` folder. After making changes to the file, you need to restart and rerun this notebook.
<a id='poissonp1.hh'> </a>


```c++
#include "poissonp1.hh"

using LOP = PoissonP1<decltype(f),FEM>;
LOP lop(f, fem.find(*gv.template begin<0>()));
```

*The following solutions are again for exercise 2*
<a id='solex2LOP'> </a>

*The cell beneath can be either executed with your own modification of the poisson class `mypoisson1.hh` or with the solution `solutionpoissonp1.hh`.*


```c++
/*#include "mypoissonp1.hh"
//#include "solutionpoissonp1.hh"

using LOP = PoissonP1<decltype(f),decltype(k),decltype(a),FEM>;
LOP lop(f,k,a,fem.find(*gv.template begin<0>()));*/
```

This local operator is used as one of the template arguments in the global assembler or grid operator, which additionally requires the types for ansatz and test space, a matrix backend, the types for the components of the coefficient vectors of ansatz and test space and the entries of the matrix *A* and finally the types of the constraints container for test and ansatz space.



```c++
using MBE = Dune::PDELab::ISTL::BCRSMatrixBackend<>;
MBE mbe(1<<(dim+1)); // guess nonzeros per row
using GO = Dune::PDELab::GridOperator<
  GFS,GFS,  /* ansatz and test space */
  LOP,      /* local operator */
  MBE,      /* matrix backend */
  RF,RF,RF, /* domain, range, jacobian field type*/
  CC,CC     /* constraints for ansatz and test space */
  >;
GO go(gfs,cc,gfs,cc,lop,mbe);
```

Next, a linear solver for the linear system $Az = b$ is selected. Since ISTL backends for vectors and matrices are used, also a linear solver from the ISTL library is chosen.
Complete solvers are packaged by PDELab for sequential and
parallel computations. Here we select the conjugate gradient method with algebraic
multigrid as preconditioner and symmetric successive overrelaxation as smoother in
multigrid in its sequential implementation. The linear solver object `ls` is initialized
with the maximum number of iterations and a verbose parameter:


```c++
using LS = Dune::PDELab::ISTLBackend_SEQ_CG_AMG_SSOR<GO>;
LS ls(100,2);
```

*The following cell contains a different linear iterative solver (see exercise 1.3). Again, only one of the cells should be executed.*


```c++
 //using LS = Dune::PDELab::ISTLBackend_SEQ_BCGS_SSOR;
 //LS ls(5000,true);
```

So far, no finite element computations have actually been performed, only the
necessary components have been configured to now do the real work together. We
can now set up the matrix *A* as well as the right hand side *b* and then solve the
system. Since this is required often, there is a class in PDELab for this:


```c++
using SLP = Dune::PDELab::StationaryLinearProblemSolver<GO,LS,Z>;
SLP slp(go,ls,z,1e-10);
slp.apply(); // here all the work is done!
```

In the given example a problem with a known exact solution, which is given by
the function *g* is solved. In order to compare the computed solution with the exact
solution we initialize another coeffcient vector with the Lagrange interpolant of the
exact solution and provide a grid function:


```c++
using ZDGF = Dune::PDELab::DiscreteGridFunction<GFS,Z>;
ZDGF zdgf(gfs,z);
Z w(gfs); // Lagrange interpolation of exact solution
Dune::PDELab::interpolate(g,gfs,w);
ZDGF wdgf(gfs,w);
Dune::VTKWriter<GV> vtkwriter(gv,Dune::VTK::conforming);
using VTKF = Dune::PDELab::VTKGridFunctionAdapter<ZDGF>;
vtkwriter.addVertexData(std::shared_ptr<VTKF>(new VTKF(zdgf,"fesol")));
vtkwriter.addVertexData(std::shared_ptr<VTKF>(new VTKF(wdgf,"exact")));
```

The visualization data can be generated in Jupyter by printing the VTKWriter instance:


```c++
vtkwriter
```

# Local Operator 

The finite element method itself is implemented in the so-called *local operator* realized by the class template
`PoissonP1`. It provides all the necessary element-local computations as described **in Subsection 3.4** and is declared as follows:

```c++
template<typename F, typename FiniteElementMap>
class PoissonP1 :

  public Dune::PDELab::FullVolumePattern,
  public Dune::PDELab::LocalOperatorDefaultFlags
```
The first template parameter provides the right hand side function of the PDE and
the second parameter provides a finite element map giving access to finite element
basis functions on the reference element for all elements of the grid. The class derives
from the PDELab classes `FullVolumePattern` and `LocalOperatorDefaultFlags` which provide some default constants and methods.

The basic assumption of this implementation of the finite element method is that all elements of the mesh are simplices of dimension d which use the same polynomial degree 1. In order to make the code faster it is a good idea to do the evaluation of the basis functions and their gradients on the reference element once before the
computations start. This will be done in the constructor, but before we can do so we need to do some preparations.

## Type Definitions and Data Members

The class begins by extracting important types. The finite element map provides a finite element for each element of the map. Its type is
```c++
private:
  // define useful types
  typedef typename FiniteElementMap::Traits::FiniteElementType
     FiniteElementType;
```     
Among other things the finite element contains the basis functions on the reference element which can be accessed via the following type:
```c++
  typedef typename FiniteElementType::Traits::LocalBasisType
     LocalBasisType;
```     
DUNE thinks of basis functions on the reference element to be of
the most general form 
$$\hat\phi : \mathbb{A}^d \to \mathbb{B}^k, \quad 
\nabla\hat\phi : \mathbb{A}^d \to \mathbb{B}^{k\times d},$$
i.e. they may be vector-valued. The following type definitions     
```c++
  typedef typename LocalBasisType::Traits::DomainType
     DomainType;
  typedef typename LocalBasisType::Traits::RangeFieldType
     RF;
  typedef typename LocalBasisType::Traits::RangeType
     RangeType;
  typedef typename LocalBasisType::Traits::JacobianType
     JacobianType;
```     
provide types to represent arguments and results of basis function
evaluations. `DomainType` represents $\mathbb{A}^d$,
`RF` represents $\mathbb{B}$, `RangeType`
represents $\mathbb{B}^k$ and finally `JacobianType`
represents $\mathbb{B}^{k\times d}$.   

Next, we extract some important constants, the dimension 
of the grid and the number of basis functions per element:
```c++
  // data members
  enum {dim=LocalBasisType::Traits::dimDomain};
  enum {n=dim+1};
```  
As private data members the class stores an instance of the right hand side function `f` provided by the `driver`:
```c++
  const F f;              // right hand side function
```  
the midpoint quadrature rule
```c++
  DomainType qp;          // center of mass of refelem
  double weight;          // quadrature weight on refelem
```  
where `qp` is $\hat S_d$ and `weight` is $w_d$, and the values of the basis functions at the quadrature
point and their gradients:
```c++
  double phihat[n];       // basis functions at qp
  double gradhat[dim][n]; // coordinate x #basisfct
```  
Then, already in the public part, we need to define some constants that control the operation of the grid operator doing the global assembly:
```c++
public:
  // define flags controlling global assembler
  enum { doPatternVolume = true };
  enum { doAlphaVolume = true };
  enum { doLambdaVolume = true };
```  
These constants are evaluated at *compile time* and tell he grid operator class which methods have been implemented in the local operator by the user. Actually, the base class `LocalOperatorDefaultFlags` provides all possible flags ith the value`false` and we just need to overwrite the nes that are needed.
The constant`doPatternVolume` ells the global assembler to determine the sparsity pattern of the
matrix $A$ from a method `pattern_volume` which i inherited from the base class `FullVolumePattern`.
This default implementation inserts nonzeros between all degrees f freedom of an element. The constants `doAlphaVolume` and `doLambdaVolume` determine that our finite element method contains a volume integral involving the finite element solution $u_h$ and a right hand side integral which does not involve the finite element solution.

Setting `doAlphaVolume` to true implies that the local operator class implements the methods `alpha_volume`, 
`jacobian_apply_volume` and `jacobian_volume`.
Setting `doLambdaVolume` to true implies that the method `lambda_volume` must be implemented.

## Constructor

The constructor of the class has the following signature:
```c++
 // Constructor precomputes element independent data
  PoissonP1 (const F& f_, const FiniteElementType& fel)
    : f(f_)
```

It takes the right hand side function `f` and a finite element `fel` as argument. The finite element is obtained from the finite element map and the first element of the grid in the function `driver`.
First thing to do is to get the lowest order quadrature rule for simplices from DUNE and we check that this is actually the midpoint rule:
```c++
  {
    // select quadrature rule
    Dune::GeometryType gt = fel.type();
    const Dune::QuadratureRule<RF,dim>&
      rule = Dune::QuadratureRules<RF,dim>::rule(gt,1);
    if (rule.size()>1) {
      std::cout << "Wrong quadrature rule!" << std::endl;
      exit(1);
    }
```
Then we store the first quadrature point in the local data members:
```c++
    // position and weight of the quadrature point
    weight = rule[0].weight();
    qp = rule[0].position();
```    
It is also a good idea to check that the basis given by the user has at least the correct size:
```c++
    // check size of the basis
    if (fel.localBasis().size()!=n) {
      std::cout << "Wrong basis!" << std::endl;
      exit(1);
    }
```    
Now the basis functions can be evaluated at the quadrature point in the reference element and the results are stored in the data members of the class:
```c++
    // evaluate basis functions on refelem
    std::vector<RangeType> phi(n);
    fel.localBasis().evaluateFunction(qp,phi);
    for (int i=0; i<n; i++) phihat[i] = phi[i];
```    
And the same now for the gradients:
```c++
    // evaluate gradients of basis functions on refelem
    std::vector<JacobianType> js(n);
    fel.localBasis().evaluateJacobian(qp,js);
    for (int i=0; i<n; i++)
      for (int j=0; j<dim; j++)
        gradhat[j][i] = js[i][0][j];
  }
```  
Note that the last index loops over the number of basis functions. 

## Method `lambda_volume`

This method computes the contributions $b_T$ to the right hand side vector for a given element. It has the following signature:
```c++
  // volume integral depending only on test functions
  template<typename EG, typename LFSV, typename R>
  void lambda_volume (const EG& eg, const LFSV& lfsv,
                      R& r) const
```                      
Argument`eg` provides the element $T$ in a wrapped form such that PDELab need not operate directly on a DUNE grid.
With `eg.geometry()` the geometry of the element can be accessed in the form of a `Dune::Geometry`. With `eg.entity()` one can access the underlying codim 0 entity of the DUNE grid.
The second argument `lfsv` provides the test functions on the reference element and `r` provides a container where the result should be stored.
First thing to do is to evaluate the right hand side function at the quadrature point:
```c++
  {
    typename F::Traits::RangeType fval;
    f.evaluate(eg.entity(),qp,fval);
```    
Next, we compute the factor that is common to all entries of $b_T$: 
```c++ 
    RF factor=fval*weight*eg.geometry().integrationElement(qp);
```  

Note that the method `integrationElement` on the geometry provides the value of $|\det B_T|$.
Finally, we can compute the entries and store them in the results container:

```c++
    for (int i=0; i<n; i++)
      r.accumulate(lfsv,i,-factor*phihat[i]);
  }
```  
Here it is important to note the minus sign because PDELab actually solves the weak formulation as 
$$r(u_h,v) = a(u_h,v)-l(v) = 0 \quad \forall v\in V$$ since this is more appropriate in the case of nonlinear partial differential equations.

## Method `jacobian_volume`

Next we need to compute the element contributions to the stiffness matrix.
This is done by the method `jacobian_volume` with
```c++
  // jacobian of volume term
  template<typename EG, typename LFSU, typename X,
           typename LFSV, typename M>
  void jacobian_volume (const EG& eg, const LFSU& lfsu,
                        const X& x, const LFSV& lfsv,
                        M& mat) const
  {
```  
Its arguments are: `eg` providing the wrapped codim 0 entity $T$,`lfsu` providing the basis functions of the ansatz space, `x` providing the coefficients of the current iterate of the finite element solution, `lfsv` providing the test functions and `mat` a container to store the result.

The `jacobian_volume` method works in the same way also for nonlinear problems. Nonlinear problems are solved iteratively, e.g. using Newton's method or a fixed-point iteration, where the method should provide a linearization at the current iterate given by the combination of `lfsu` and `x`.
In our case of a linear problem the result *does not depend* on the current iterate.
Moreover, the basis functions for the test space are precomputed so we need
not access them via `lfsv`. Note also that in general the ansatz and test space might be different.

First thing we need is to get $B_T^{-T}$ and store it into `S`:
```c++
    // get Jacobian and determinant
    // assume the transformation is linear
    const auto geo = eg.geometry();
    const auto S = geo.jacobianInverseTransposed(qp);
```    
Next, $|\det B_T|$ is retrieved from the geometry and the factor that is common to all entries of the local stiffness matrix is computed:
```c++
RF factor = weight*geo.integrationElement(qp);
```
Now form the matrix of transformed gradients $G=B_T^{-T} \hat G$ and store it in `grad`:
```c++
  // compute gradients of basis functions in transformed element
    double grad[dim][n] = {{0.0}}; // coordinate x #basisfct
    for (int i=0; i<dim; i++) // rows of S
      for (int k=0; k<dim; k++) // columns of S
        for (int j=0; j<n; j++) // columns of gradhat
          grad[i][j] += S[i][k] * gradhat[k][j];
```          
The computations are arranged in such a way that the innermost loop has the dimension number of basis functions. In 3d there are four basis functions and the loop has a chance to get vectorized.
 
Now the local stiffness matrix $A_T = G^T G$ (up to the factor $|\det B_T| w_d$) is formed
```c++
    // compute grad^T * grad
    double A[n][n] = {{0.0}};
    for (int i=0; i<n; i++)
      for (int k=0; k<dim; k++)
        for (int j=0; j<n; j++)
          A[i][j] += grad[k][i]*grad[k][j];
```          
and stored in the results container (now multiplying with the common factor):
```c++
    // store in result
    for (int i=0; i<n; i++)
      for (int j=0; j<n; j++)
        mat.accumulate(lfsu,i,lfsu,j,A[i][j]*factor);
  }
```

## Method `alpha_volume`

The method `alpha_volume` provides the element-local computations for the matrix-free evaluation of $a(u_h,\phi_i)$ for all test functions $\phi_i$. It has the interface:
```c++
 // volume integral depending on test and ansatz functions
  template<typename EG, typename LFSU, typename X,
           typename LFSV, typename R>
  void alpha_volume (const EG& eg, const LFSU& lfsu,
                     const X& x, const LFSV& lfsv,
                     R& r) const
```                     
Its arguments are: `eg` providing the wrapped codim 0 entity $T$,`lfsu` providing the basis functions of the ansatz space,`x` providing the coefficients of the current iterate of the finite element solution, `lfsv` providing the test functions and `r` a container to store the result. 

The computations are actually quite similar to those in `jacobian_volume`. In particular, the computation of $B_T^{-T}$, $|\det B_T|$ and $G = B_T^{-T} \hat G$ are the same.
```c++
  {
    // get Jacobian and determinant
    // assume the transformation is linear
    const auto geo = eg.geometry();
    const auto S = geo.jacobianInverseTransposed(qp);
    RF factor = weight*geo.integrationElement(qp);

    // compute gradients of basis functions in transformed element
    double grad[dim][n] = {{0.0}};  // coordinate x #basisfct
    for (int i=0; i<dim; i++) // rows of S
      for (int k=0; k<dim; k++) // columns of S
        for (int j=0; j<n; j++) // columns of gradhat
          grad[i][j] += S[i][k] * gradhat[k][j];
```          
Extracting the element local coefficients $z_T = R_T z$ is done by:
```c++
    double z_T[n];
    for (int j=0; j<n; j++) z_T[j] = x(lfsu,j); // read coeffs
```    
Now we may compute $\nabla u_h$ via $G z_T$:
```c++
    double graduh[dim] = {0.0};
    for (int k=0; k<dim; k++) // rows of grad
      for (int j=0; j<n; j++) // columns of grad
        graduh[k] += grad[k][j]*z_T[j];
```        
Finally, the result $a_T = G^T \nabla u_h$ is formed:
```c++
    // scalar products
    double a_T[n] = {0.0};
    for (int k=0; k<dim; k++) // rows of grad
      for (int j=0; j<n; j++)
        a_T[j] += grad[k][j]*graduh[k];
```
and stored in the results container (while being multiplied with the common factor):
```c++
    // store in result
    for (int i=0; i<n; i++)
      r.accumulate(lfsv,i,a_T[i]*factor);
  }
```

## Method `jacobian_apply_volume`

In the case of a *nonlinear* partial differential equation the finite element method results in a weak form 
$$u_h \in V_h \ : \quad r(u_h,v) = \alpha(u_h,v) - \lambda(v) = 0 \quad \forall v\in V_h$$
which is *nonlinear* in its *first* argument. Inserting the finite element basis results in a nonlinear algebraic problem
$$R(z) = 0$$
with $(R(z))_i = r(\text{FE}_h(z),\phi_i)$ which is typically solved by Newton's iteration
or some other iterative method. In case of Newton's method, each step
involves the solution of a {\em linear} system of the form
$$J(z)\, w = R(z)$$
where $(J(z))_{i,j} = \frac{(\partial R(z))_i}{\partial z_j} = 
\frac{\partial \alpha(\text{FE}_h(z),\phi_i)}{\partial z_j}$ is the Jacobian of 
the nonlinear map $R$.

Naturally, the nonlinear case also includes the linear case described in this tutorial
by setting $r(u,v) = a(u,v)-l(v)$. Then, due to the linearity of $a$ in its first
argument, one can show that $J(z) = A$ and 
$$(J(z)\, w)_i = (Aw)_i = a(\text{FE}_h(w),\phi_i).$$
This is *not* true in the nonlinear case. There, the evaluation of the form $\alpha(\text{FE}_h(w),\phi_i)$ and the application of the Jacobian $J w$ are different operations. Therefore, PDELab provides two functions with
the application of the Jacobian implemented in `jacobian_apply_volume`
with the interface:
```c++
  //! apply local jacobian of the volume term
  template<typename EG, typename LFSU, typename X,
           typename LFSV, typename R>
  void jacobian_apply_volume (const EG& eg, const LFSU& lfsu,
                              const X& z, const LFSV& lfsv,
                              R& r) const
```                              
Note, this is the same interface as for `alpha_volume`. Since our problem is linear, Jacobian application is identical to bilinear form evaluation and therefore we may just forward the call to the
function `alpha_volume`: 

```c++
  {
    alpha_volume(eg,lfsu,z,lfsv,r);
  }
``` 

# Exercises
<a id='Exercises'> </a>

## Exercise 1

1. **Warming up**

   Run the program with different refinement levels. Visualize the solution in ParaView, use the Calculator
   filter to visualize the difference $|u-u_h|$ and determine the
   maximum error. Note that $u$ and $u_h$ are called `exact` and `fesol` in the paraview output.

- *In the following $|u-u_h|$ for different refinement levels is depicted.*

<table  style="width:100%">
  <tr>
    <td></td>
    <td></td>
  </tr>
  <tr>
    <td></td>
    <td></td>
  </tr>

2. **Solving a new problem**

   Now consider
   \begin{equation}
   f(x) = -\sum_{i = 1}^{d} 6(x)_i \quad \text{and} \quad g(x) = \sum_{i = 1}^d(x)_i^3
   \end{equation}
   and check that $u(x)=\sum_{i=1}^d (x)_i^3$ solves the PDE. Implement the new $f$ and $g$ and rerun your program.

- *The solution for the new implementation of f and g is embedded in the notebook, directly at the place of the defintions ([go to solution](#newfandg)). When running all cells, ensure that only one of the definitions is uncommented (a cell that is hidden through the 'Show/Hide Solution' button is still executed).*

  

3. **Analysis of finite element error**:

  Produce a sequence of output files for different levels of mesh
  refinements $0, 1, 2, \ldots$ with suitable output
  filenames. Visualize the error $|u-u_h|$ in Para\-View and determine
  the maximum error on each level.  If you want you can also try to
  calculate the L2-norm of the error in paraview and determine the
  convergence rate by hand. In a later exercise you will see how to
  calculate the L2-norm directly in the c++ source.

 - *The following table shows the maximum norm,which is given as $||f||_{\max} := \max|f(x)|$ and the L2-norm as $||f||_{0,\Omega} = \left( \int_{\Omega} f(x)^2 dx \right)^{1/2} $ for different mesh refinements. The L2 error is depicted in the graph on the right side.* 

<table  style="width:100% ">
<tr>
<td>
<table  style="border:1px dotted; width:100%">
  <tr style="border:1px dotted">
    <th>refinement</th>
    <th>${||u-u_h||}_\max$</th>
    <th>${||u-u_h||}_{0,\Omega}$</th>
  </tr>
  <tr>
    <td>0</td>
    <td>1.18978e-05</td>
    <td>1.161176e-03</td>
  </tr>
  <tr>
    <td>1</td>
    <td>1.76823e-06</td>
    <td>2.59592e-04</td>
  </tr>
  <tr>
    <td>2</td>
    <td>2.02272e-07</td>
    <td>6.0531e-05</td>
  </tr> 
  <tr>
    <td>3</td>
    <td>1.9762e-08</td>
    <td>1.4689e-05</td>
  </tr>
  <tr>
    <td>4</td>
    <td>1.7708e-09</td>
    <td>3.635e-06</td>
  </tr>
  <tr>
    <td>5</td>
    <td>1.50032e-10</td>
    <td>9.06e-07</td>
  </tr>     
</table> 
</td>

<td>
    
</td>
    
 </tr>

 - *In general the convergence behaviour of the finite element method is
$$ {||u-u_h||}_{j,\Omega} \rightarrow 0 \quad \text{for } h \rightarrow 0,$$
with $j = 0$ for the L2-norm. It can be profen that
$${||u-u_h||}_{j,\Omega} \leq Ch^{\beta},$$
where the exponent $\beta$ depends on the polynomial degree, the norm in which the error is measured and the regularity of the solution (more details can be found [here](#FurtherReading)). The exponent $\beta$ is called the* convergence rate *of the method (w.r.t a given norm). With $\beta = k + 1 $, where k denotes the polynomial degree, the finite element method with polynomials of degree 1 has a convergence rate of $\beta = 2$ w.r.t th L2-norm. The plot on the right depicts the error in the L2-norm with respect to the mesh size, clearly indicating a quadratic convergence rate.*

 - *These are the basic steps in order to determine the L2-norm of the error in paraview/by hand:*   
     1. *use the calculator for $(u(x)-u_h(x))^2$*
     2. *use 'integrate variables' filter for $\int_{\Omega}(u(x)-u_h(x))^2 \,dx$*
     3. *determine $\left( \int_{\Omega}(u(x)-u_h(x))^2 \right)^{1/2}$ by hand*

4. **Use a different solver**

   You may also try to exchange the iterative linear solver, the lines for the new solver are already given as solution next to the old linear solver. Compare the number of iterations which is given by the following lines
  in the output (here we have 12 iterations):


```bash
=== CGSolver

12      1.9016e-09

=== rate=0.144679, T=0.02362, TIT=0.00196833, IT=12
```

*The previous solver uses 7 iterations and the new solver uses 10 iterations.*

## Exercise 2

Now consider the extended equation of the form
\begin{align}
    \begin{array}{rcll}
      -\nabla\cdot (k(x) \nabla u) +a(x) u  & = & f & \text{ in } \Omega, \\
      u & = & g & \text{ on } \partial\Omega,
    \end{array}
\end{align}
with scalar functions $k(x)$, $a(x)$.

This is the linear convection-diffusion equation, with $k(x)$ being the diffusion coefficient and $a(x)$ the reaction term.

1. In a first step show that the weak formulation of the problem
  involves the new bilinear form
  \begin{align}
    a(u,v) = \int_\Omega k(x) \nabla u(x) \cdot \nabla v(x) + a(x) u(x) v(x) \,dx .
  \end{align}
<a id='newBilinearForm'> </a>

*This closely follows the derivation of the bilinear form as shown in section [Weak Formulation](#weakFormulation)  for the unextended equation.*

2. The main work is to extend the local operator given by the class
  `PoissonP1`. Your modifications should be done in the file `mypoissonp1.hh`, which is simply a copy of `poissonp1.hh`. The extension can be done in several steps:

 - Provide analytic functions for $k(x)$ and $a(x)$ and pass
    grid functions to the local operator like it is done for $f(x)$.
 - Extend the local operator to first handle $k(x)$ by assuming the
    function $k(x)$ to be *constant* on mesh elements.
 - Now extend the local operator (in `mypoissonp1.hh`) to handle $a(x) u(x)$. Also assume the function $a(x)$ to be *constant* on mesh elements.


- *The solution for providing the analytic functions and passing grid functions to the local operator are embedded in the notebook [(go to solution)](#solex2)*. 
- *The place to include your modified local operator or the solution can be found [here](solex2LOP)*.
- *The solution for changes within the local operator can be found in `solutionpoissonp1.hh`. Additionally, in the following an overview of the main changes within the local operator is provided*.

Overview Changes in `poissonp1.hh`
 - *The local operator is implemented as the class PoissonP1. In order to extend the local operator to handle the functions $k(x)$ and $a(x)$, it is provided with two additional data members, k and a. Note, that these are **scalars**, as $k(x)$ and $a(x)$ are constant on mesh elements. As these are provided from outside the local operator, the constructor needs to be adapted as well.*

```cpp
// data members
  ...
  const K k;              // diffusion coefficient
  const T a;              // reaction term

// Constructor precomputes element independent data
  PoissonP1 (const F& f_, const K& k_, const T& a_, const FiniteElementType& fel)
    : f(f_), k(k_), a(a_)
```

- *As described before, the method `jacobian_volume` computes the element contributions to the stiffness matrix for a given element. As the [new bilinear form $a(u,v)$](#newBilinearForm) includes $a(x)$ and $k(x)$, these need to be included in `jacobian_volume` as well.*

```cpp
void jacobian_volume (...) const 
{
    // evaluate diffusion coefficient
    typename K::Traits::RangeType kval;
    k.evaluate(eg.entity(),qp,kval);

    // evaluate reaction coefficient
    typename T::Traits::RangeType aval;
    a.evaluate(eg.entity(),qp,aval);           
               ...                
    mat.accumulate(lfsu,i,lfsu,j,(kval*A[i][j]+aval*phihat[j]*phihat[i] )*factor);
}
```

- *Hence also the method `alpha_volume` needs to be adadpted. As a reminder, it provides the element-local computations for the matrix-free evaluation of $a(u_h,\phi_i)$ for all test functions $\phi_i$.*

```cpp
void alpha_volume (...) const{
    // evaluate diffusion coefficient
    typename K::Traits::RangeType kval;
    k.evaluate(eg.entity(),qp,kval);

    // evaluate reaction coefficient
    typename T::Traits::RangeType aval;
    a.evaluate(eg.entity(),qp,aval);           
    ...
    r.accumulate(lfsv,i,(kval*a_T[i]+aval*u*phihat[i])*factor);
}
```

3. Use your new local operator and test it. In the solution the code is tested for $k(x) = (x)_1$,
  $a(x)=\sum_{i=1}^{d}(x)_i$ and $g=\sum_{i=1}^d (x)_i^3$ by choosing  the right hand side $f(x)$ in such a way that $u(x)=\sum_{i=1}^d (x)_i^3$ solves the PDE.

 - *Reversed engineering yields $f(x) = - \left[6 \left(\sum_{i=1}^{d}(x)_i - (x)_1\right)+ 9(x)_1 \right] \cdot x_1 + \sum_{i=1}^{d}(x)_i \cdot \sum_{j=1}^{d}(x)^3_i$.* 
 
*Below the solution is depicted.*



# Further Reading

Finite element method:
1. K. Eriksson, D. Estep, P. Hansbo and C. Johnson. *Computational Differential Equations*. Cambridge University Press.1996. 
2. A. Ern and J.-L. Guermond. *Theory and practice of finite element methods*. Springer. 2004.
3. P. G. Ciarlet. *The finite element method for elliptic problems*. SIAM. Classics in Applied Mathematics. 2002.
4. D. Braess. *Finite Elemente*. Springer. 2003.
5. S. C. Brenner and L. R. Scott. *The mathematical theory of finite element methods*. Springer. 1994.
6. H. Elman, D. Silvester and A. Wathen.*Finite Elements and Fast Iterative Solvers*. Oxford University Press. 2005.
7. C. Großmann and H.-G. Roos. *Numerische Behandlung partieller Differentialgleichungen*. Teubner. 2006.
8. W. Hackbusch. *Theorie und Numerik elliptischer Differentialgleichungen*. Teubner. 1986.
9. R. Rannacher. *Einführung in die Numerische Mathematik II (Numerik partieller Differentialgleichungen)*. 2006. http://numerik.iwr.uni-heidelberg.de/~lehre/notes.
10. P. Bastian. *Lecture Notes on Scientific Computing with Partial Differential Equations*. 2014. https://conan.iwr.uni-heidelberg.de/data/teaching/finiteelements_ws2017/num2.pdf.

<a id='FurtherReading'> </a>

9e96baa70008c4b464049bf05d437bc25ed4cac4
ipynb
Jupyter Notebook
notebooks/tutorial00/pdelab-tutorial00.ipynb
dokempf/dune-jupyter-course
1da9c0c2a056952a738e8c7f5aa5aa00fb59442c
BSD-3-Clause
2022-01-21T03:16:12.000Z
2022-01-21T03:16:12.000Z
notebooks/tutorial00/pdelab-tutorial00.ipynb
dokempf/dune-jupyter-course
1da9c0c2a056952a738e8c7f5aa5aa00fb59442c
BSD-3-Clause
2021-04-22T13:52:59.000Z
2021-10-04T13:31:59.000Z
notebooks/tutorial00/pdelab-tutorial00.ipynb
dokempf/dune-jupyter-course
1da9c0c2a056952a738e8c7f5aa5aa00fb59442c
BSD-3-Clause
2021-04-21T08:20:02.000Z
2021-04-21T08:20:02.000Z
Qwen/Qwen-72B
1. YES
2. YES