Specializations

Now, we have a basic FEM solver with rudimentary features. We can now start to add more features to the solver, in mini-projects that are done in smaller groups in parallel.

Each project should contain some new feature implementation, a demo notebook and corresponding unit tests. The feature should be implemented in a way that it could be used in the final (merged) project by all other users.

The following list is a colletion of ideas and partially pretty roughly sketched only:

Autodiff for Derivatives [ Johann Meyer + Samira Altpeter]

We have used numerical differentiation as a fall back to compute the derivatives of shape functions. Alternatively we can directly compute the derivatives alongside the shapes itself using an autodiff data type.

To this end we need to introduce a datatype which re-implements standard arithmetic operations. It holds a value and a derivative value. The arithmetic operations then compute value and derivative using standard rules of differentiation, e.g. product rule:

x = AD(ip[0],(1,0)) # the expression corresponding to x has value ip[0] and derivative (1,0)
y = AD(ip[1],(0,1)) # the expression corresponding to y has value ip[1] and derivative (0,1)
z = x*y # the expression corresponding to z
# z.val == x.val * y.val
# z.der == x.val * y.der + x.der * y.val

The task of this mini-project is to implement such a datatype and use it to compute the derivatives of shape functions.

Tasks:

• Implement AD datatype for "forward accumulation AD" (including tests, demo notebook)
• Use AD datatype for the implementation of shape functions (including tests, demo notebook)

Results:

The results of this project can be found in this summary Markdown file. The realization is then found in the autodiff branch.

MeshFunction arithmetics [ Ilya Gutkin + Maximilian Mogk]

Given several MeshFunctions (e.g. exact solution and discrete solution, i.e. a GlobalFunction and a FEFunction) we want to be able to compute the sum, difference, product, ... of these MeshFunctions as a new MeshFunction. This requires to implement the corresponding arithmetics for MeshFunctions.

• Overload arithmetic operations like __add__, __sub__, __mul__, __div__, __invert__, __neg__, __pos__ for MeshFunctions (including tests, demo notebook)
• Use them to compute the difference between exact and discrete solution (including tests, demo notebook)
• Provide a gradient operator for MeshFunctions (non-implemented for general MeshFunctions but implemented for FEFunctions) (including tests, demo notebook)

Results:

The results of this project can be found in this notebook. The realization is then found in the convdiff branch.

Consider a different scalar PDE: Convection-Diffusion equations [Mika Meyer, Kolja Straatman, Pascal Irmer]

Consider the (reaction-)convection-diffusion equation $$ -\nabla \cdot (\alpha \nabla u) + \beta \cdot \nabla u + \gamma u = f $$ in $\Omega$, equipped with suitable boundary conditions.

Implement a FEM solver for this problem and carry out numerical experiments to investigate the influence of the convection term on the solution.

Consider the following standard test problem:

• Set $ \gamma = 0$, and $\beta = (1,1)$, $\Omega = (0,1)^2$
• Use homogeneous Dirichlet boundary conditions: $u|_{\partial \Omega} = 0$, see Dirichlet demo notebook.
• Consider the so-called mesh Péclét number $Pe = \frac{2 \| \beta \|_\infty}{\gamma h}$, where $h$ is the (initially uniform) mesh size.
• Take the details as in this paper, page 19 ff., i.e. take the exact solution and compute (verify) the corresponding r.h.s. and compute the numerical solutions.
• Implement the (missing) convection integral for $\int_\Omega \beta \cdot \nabla u v \, dx $ in the bilinear form (including tests and notebook(s))
• Do numerical studies over $h$ and/or $Pe$ with different FE spaces
• Document the results in proper notebook(s)
• Use strechted grids (by introducing a mapping in the construction of the mesh) with increasingly small elements towards $x=1$ and $y=1$ and investigate the influence on the solution. (including documentation in notebook(s))

Results:

The results of this project can be found in this notebook. The realization is then found in the adaptive branch.

Adaptive mesh refinement [ Paula John + Marvin Langer ]

Implement an adaptive mesh refinement strategy for the Poisson problem. This requires to implement a posteriori error estimation and a mesh refinement strategy.

Tasks:

• Consider a 1D problem of the form $- \alpha u'' + \gamma u = f$ with $\gamma = 1$, $\alpha \ll 1$ and $f = \chi_{[0.4,0.6]}$ (i.e. a bump in the middle of the domain), $\Omega = (0,1)$.

• Compute the numerical solution for $\alpha \in \{ 10^{-i} \}_{i=1,2,3,4,5}$. What do you observe?

• After computation of the discrete solution, implement an a posteriori error estimation based on the following error indicator that is evaluated on each element $T \in \mathcal{T}_h$:

  $$ \eta_T = h_T^2 \| - \alpha u_h'' + \gamma u_h - f \|_{L^2(T)}^2 + \sum_{F \subset \partial T} h_T \| \llbracket u_h' \rrbracket \|_{L^2(F)}^2 $$

  Here $h_T$ and $h_F$ denote the diameter of the element $T$ and the face $F$, respectively, and $\llbracket u' \rrbracket$ denotes the jump of the derivative of $u$ across the face $F$ (vertex in 1D case).

• Visualize the error indicator on the mesh (as a piecewise constant function).

• Use a "Dörfler" marking, i.e. select the elements with largest $\eta_T$ that add up to $20 \%$ of the total error indicator $\sum_T \eta_T$.

• Visualize the marked mesh.

• Now refine the mesh according to the marking and recompute the solution. Repeat this until the (summed) error indicator is below a certain tolerance.

• Document and discuss the new implementation and the results in a notebook.

• Think about extensions to 2D. Sketch a possible implementation based on a longest edge or a red-green refinement strategy. Explain the additional challenges that arise in 2D.

Results:

The results of this project can be found in this notebook. The realization is then found in the adaptive branch.

Consider different linear solvers [ Maximilian Zienecker]

So far we only used a direct solver to solve the PDE (which is not that bad in 1D and 2D, but not optimal in 3D). In this task we want to consider different linear solvers, e.g. a CG solver with preconditioning.

• Consider the Poisson problem in 2D (with Dirichlet boundary conditions or a mass term)
• Use a CG solver without preconditioning to solve the linear systems (including tests, demo notebook). In the notebooks also explain the numpy/scipy functionality used.
• Investigate standard scipy-preconditioners. Compare number of iterations and runtime (including tests, demo notebook)
• Implement a simple preconditioners and compare to the standard preconditioners (including tests, demo notebook). Simple preconditioner choices are:
  • Jacobi preconditioner (diagonal of the matrix)
  • Block-Jacobi preconditioner (diagonal blocks of the matrix, with blocks according to a partition of the mesh)
  • Overlapping block Jacobi preconditioner (diagonal blocks of the matrix, with blocks according to a partition of the mesh, but with overlap, e.g. all elements that share a vertex in one block)

Results:

The results of this project can be found in this notebook. The realization is then found in the linsolve branch.

---

The following projects have not been touched within the course so far:

Meshes and Finite elements on quadrilaterals

Extend mesh structure and finite elements to quadrilaterals. This includes the following tasks:

• Implement a simple 2D mesh based on quads
• Implement quadrilateral finite elements (based on tensor products of 1D elements)
• Implement finite element spaces (or extend existing ones) for quadrilateral

Higher order (non-Lagrange) FE functions in 2D

Higher order basis functions can be constructed based on tensor product finite elements for the quadrilateral in combination with a Duffy transformation, resulting in 2D polynomials based on Jacobi and integrated Legendre polynomials.

Boundary integrals

So far, all considered integrals have been integrals on $\Omega$. However, other use cases are integrals on $\partial \Omega$. For example, the boundary integral that arises from non-homogeneous Neumann boundary conditions.

The task of this mini-project is to implement a boundary integral:

• This requires to change the assembly-functions to distinguish between loops over boundary and volume elements and to apply corresponding FormIntegrators only the corresponding parts
• This again requires to change the FormIntegrator interface to allow for different integrals on boundary and volume elements.
• Further a transformation from 1D to 2D, respectively 0D to 1D is required to map the boundary integral to the reference element (of co-dimension 1).

Non-homogeneous Dirichlet boundary conditions

So far, we have only implemented homogenenous Dirichlet boundary conditions (and potentially only with a hack). In this task we want to implement non-homogeneous Dirichlet boundary conditions more systematically in several steps:

• Projection of a given function on a given set of boundary elements, using the Oswald interpolation operator which applys a simple $L^2$ projection on all relevant boundary elements. For dofs that appear several times the corresponding mean value is used. Hence, here, the $L^2$ projection needs to be implemented.

• Given the projection, a particular solution for a homogenization strategy can be computed. This is then used to setup a linear system for a correction PDE solution that satisfies homogeneous Dirichlet boundary conditions.

For this project further details will be given.

Oswald interpolation

Implementation of an Oswald projection operator:

• Apply a simple $L^2$ projection on all relevant boundary elements.
• For dofs that appear several times the corresponding mean value is used.

Hence, here, the $L^2$ projection needs to be implemented.

Nodal interpolation operator

Implementation of a nodal interpolation operator for Lagrange-type FESpaces.