Note

Click here to download the full example code

This tutorial covers the first steps into finite element computation
using the *M* (Modelling) in *GIMLi*.

We will not dig into deep details about the theory of finite elements
here, as this can be found in several books, e.g.,
:cite:`Zienkiewicz1977`

.

However, there is a little need for theory to understand what it means to use the finite elements method to solve a modelling problem, so we start with some basics.

Assuming the Poisson equation for a simple partial differention equation that needs to be solved for the unknown scalar field \(u(\mathbf{r})\) within a modelling domain \(\mathbf{r}\in\Omega\) for a non zero right hand side function \(f\).

\[\begin{split}- \Delta u & = f \quad{\mathrm{in}}\quad~\Omega\\
u & = g \quad{\mathrm{on}}\quad\partial\Omega\end{split}\]

The Laplace operator \(\Delta = \nabla\cdot\nabla\) given by the divergence of the gradient is the sum of second partial derivatives the of the field \(u(\mathbf{r})\) with respect to the Cartesian coordinates in 1d space \(\mathbf{r} = (x)\), in 2d \(\mathbf{r} = (x, y)\), or 3d space \(\mathbf{r} = (x, y, z)\). On the boundary \(\partial\Omega\) of the domain we want known values of \(u=g\) as the so called Dirichlet boundary conditions.

An approximated solution \(u_h\approx u\) for our numerical problem will probably only satisfy with a rest \(R\): \(\Delta u_h + f = R\). If we choose some weighting functions \(w\), we can try to minimize these residuum over our modelling domain.

\[\begin{split}\int_{\Omega} R w &= 0 \\
\int_{\Omega} - \Delta u w & = \int_{\Omega} f w\end{split}\]

It is preferable to eliminate the second derivative in the Laplace operator either due to integration by parts or by applying production rule and Gauss’s law. This leads to the so called weak formulation:

\[\begin{split}\int_{\Omega} \nabla u \nabla w - \int_{\partial \Omega}\mathbf{n}\nabla u w & = \int_{\Omega} f w \\
\int_{\Omega} \nabla u \nabla w & = \int_{\Omega} f w + \int_{\partial \Omega}\frac{\partial u}{\partial\mathbf{n}} w\end{split}\]

to be continued soon …

```
import pygimli as pg
import numpy as np
import matplotlib.pyplot as plt
```

We create a grid for our modelling domain with equidistant spacing in x and y direction.

```
grid = pg.createGrid(x=np.linspace(-1.0, 1.0, 10),
y=np.linspace(-1.0, 1.0, 10))
```

Now we can call the solver `pygimli.solver.solve`

for some
default material values and global homogeneous Dirichlet boundary conditions.

```
u = pg.solver.solve(grid, f=1.,
bc={'Dirichlet': [grid.findBoundaryByMarker(1, 5), 0.0]},
verbose=True)
```

Out:

```
Mesh: Mesh: Nodes: 100 Cells: 81 Boundaries: 180
('Asssemblation time: ', 0.002)
('Solving time: ', 0.001)
```

The result is drawn with the function `pygimli.show`

.

```
ax, cbar = pg.show(grid, data=u, label='P1 Solution $u$')
```

`pygimli.show`

is just a shortcut for various routines that can also
be called directly.

```
pg.mplviewer.drawMesh(ax, grid)
```

We repeat the computation with a spatially (H) refined version of the original grid.

```
gridh2 = grid.createH2()
uh = pg.solver.solve(gridh2, f=1.,
bc={'Dirichlet': [gridh2.findBoundaryByMarker(1, 5), 0.0]},
verbose=True)
ax, cbar = pg.show(gridh2, data=uh, label='H2 Solution $u$')
pg.mplviewer.drawMesh(ax, gridh2)
```

Out:

```
Mesh: Mesh: Nodes: 361 Cells: 324 Boundaries: 684
('Asssemblation time: ', 0.009)
('Solving time: ', 0.002)
```

We do the same using quadratic (P) refinement, i.e. the same number of nodes.

```
gridp2 = grid.createP2()
up = pg.solver.solve(gridp2, f=1.,
bc={'Dirichlet': [gridp2.findBoundaryByMarker(1, 5), 0.0]},
verbose=True)
```

Out:

```
Mesh: Mesh: Nodes: 280 Cells: 81 Boundaries: 180
('Asssemblation time: ', 0.004)
('Solving time: ', 0.002)
```

Fortunately there exist an analytical solution for this example.

To compare the different results the in detail we interpolate our solution along a probe line through the domain.

```
x = np.linspace(-1.0, 1.0, 100)
probe = np.zeros((len(x), 3))
probe[:, 0] = x
uH1 = pg.interpolate(srcMesh=grid, inVec=u, destPos=probe)
uH2 = pg.interpolate(srcMesh=gridh2, inVec=uh, destPos=probe)
uP2 = pg.interpolate(srcMesh=gridp2, inVec=up, destPos=probe)
plt.figure()
plt.plot(x, np.array(list(map(uAna, probe))), 'black', linewidth=2,
label='analytical')
plt.plot(x, uH1, label='linear (H1)')
plt.plot(x, uH2, label='linear (H2)')
plt.plot(x, uP2, label='quadratic (P2)')
plt.xlim([-0.4, 0.4])
plt.ylim([0.25, 0.3])
plt.legend()
plt.show()
```