# Introduction to FEM Analysis with Python¶

This tutorial aims to show using Python to pre-processing, solve, and post-processing of Finite Element Method analysis. It uses a finite element method library with a Python interface called GetFEM for preprocessing and solving. We will load vtk file by using meshio and visualize by matplotlib in pre-processing and post-processing. This tutorial was used in the PyConJP 2019 talk. You can watch the talk on YouTube below. This tutorial is based on the following official GetFEM page tutorial.

```
[1]:
```

```
from IPython.display import YouTubeVideo
YouTubeVideo("6JuB1GiDLQQ", start=512)
```

```
[1]:
```

## Installation¶

GetFEM including its python interface can be installed from a terminal by executing aptitude update and aptitude install python3-getfem++.

sudo aptitude install python3-getfem++The additional packages in requirements.txt are required for this tutorial. You do not need to build these environments because they are already configured in the Dockerfile.

## The problem setting¶

The problem refers to “Poisson’s Equation on Unit Disk” published by Math Works’s homepage.

## How to use GetEM¶

We take the following steps when using GetFEM to solve finite element problems. See this page for more information on using GetFEM.

- define a MesherObject
- define a Mesh
- define a MeshFem
- define a MeshIm
- define a Model and set it up
- solve Model object
- get value from Model object

## Initialization¶

GetFEM can be imported as follows (numpy has also to be imported).

```
[2]:
```

```
import getfem as gf
import numpy as np
```

## Mesh generation¶

We use GetFEM’s MesherObject to create a mesh from the geometric information to be analyzed. This object represent a geometric object to be meshed by the experimental meshing procedure of GetFEM. We can represent a disk by specifying a radius, a 2D center, and using “ball” geometry.

```
[3]:
```

```
center = [0.0, 0.0]
radius = 1.0
mo = gf.MesherObject("ball", center, radius)
```

We can make mesh object `mesh`

by calling the experimental mesher of GetFEM on the geometry represented by `mo`

. The approximate element diameter is given by `h`

and the degree of the mesh is `K`

(\(K > 1\) for curved boundaries).

```
[4]:
```

```
h = 0.1
K = 2
mesh = gf.Mesh("generate", mo, h, K)
```

## Boundary selection¶

To define a boundary condition, we set a boundary number on the outer circumference of the circle.

```
[5]:
```

```
outer_faces = mesh.outer_faces()
OUTER_BOUND = 1
mesh.set_region(OUTER_BOUND, outer_faces)
```

## Mesh draw¶

We visualize the created mesh to check its quality. We can output mesh objects, but matplotlib can only display triangles. Therefore, we make a `Slice`

object with a slice operation of `("none",)`

, which does not cut the mesh. With a refinement of 1, this serves to convert the mesh to triangles.

```
[6]:
```

```
sl = gf.Slice(("none",), mesh, 1)
```

We can export a slice to VTK file by using export_to_vtk method.

```
[7]:
```

```
sl.export_to_vtk("sl.vtk", "ascii")
```

We can render VTK files using Paraview or mayavi2. In order to display in the jupyter notebook this time, we read in meshio and draw in matplotlib.

```
[8]:
```

```
import pyvista as pv
from pyvirtualdisplay import Display
display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("sl.vtk")
p.add_mesh(m, show_edges=True)
pts = m.points
p.show(window_size=[512, 384], cpos="xy")
display.stop()
```

```
[8]:
```

```
<pyvirtualdisplay.display.Display at 0x7f96ba6ef160>
```

## Definition of finite element methods and integration method¶

We define the finite element and integration methods to use. We create a `MeshFem`

that defines the degree of freedom of the mesh as rank 1 (scalar).

```
[9]:
```

```
mfu = gf.MeshFem(mesh, 1)
```

Next we set the finite element used. Classical finite element means a continuous Lagrange element. Setting `elements_degree`

to 2 means that we will use quadratic (isoparametric) elements.

```
[10]:
```

```
elements_degree = 2
mfu.set_classical_fem(elements_degree)
```

The last thing to define is an integration method `mim`

. There is no default integration method in GetFEM so it is mandatory to define an integration method. Of course, the order of the integration method has to be chosen sufficient to make a convenient integration of the selected finite element method. Here, the square of `elements_degree`

is sufficient.

```
[11]:
```

```
mim = gf.MeshIm(mesh, pow(elements_degree, 2))
```

## Model definition¶

The model object in GetFEM gathers the variables of the models (the unknowns), the data and what are called the model bricks. The model bricks are some parts of the model (linear or nonlinear terms) applied on a single variable or linking several variables. A model brick is an object that is supposed to represent a part of a model. It aims to represent some integral terms in a weak formulation of a PDE model. They are used to make the assembly of the (tangent) linear system (see The model object for more details).

```
[12]:
```

```
md = gf.Model("real")
md.add_fem_variable("u", mfu)
```

## Poisson’s equation¶

To define Poisson’s equation, we have to define a Laplacian term and RHS source term. We can add the Laplacian term (which is called a brick in GetFEM) by using `add_Laplacian_brick()`

.

```
[13]:
```

```
md.add_Laplacian_brick(mim, "u")
```

```
[13]:
```

```
0
```

If you want to define constants in GetFEM, we use the `add_fem_data()`

method.

```
[14]:
```

```
F = 1.0
md.add_fem_data("F", mfu)
```

We can set constant values with the `set_variable()`

method. Here we pass a vector (`ndarray`

) the size of the degrees of freedom.

```
[15]:
```

```
md.set_variable("F", np.repeat(F, mfu.nbdof()))
```

We define the term RHS with the `add_source_term_brick()`

method, using the constant `F`

just defined.

```
[16]:
```

```
md.add_source_term_brick(mim, "u", "F")
```

```
[16]:
```

```
1
```

Finally, we set the Dirichlet condition at the boundary.

```
[17]:
```

```
md.add_Dirichlet_condition_with_multipliers(mim, "u", elements_degree - 1, OUTER_BOUND)
```

```
[17]:
```

```
2
```

## Model solve¶

Once the model is correctly defined, we can simply solve it by:

```
[18]:
```

```
md.solve()
```

```
[18]:
```

```
(0, 1)
```

## Export/visualization of the solution¶

The finite element problem is now solved. We can get the solution \(u\) by using `variable`

method.

```
[19]:
```

```
U = md.variable("u")
```

We can output the computed u with the mesh of the Slice object.

```
[20]:
```

```
sl.export_to_vtk("u.vtk", "ascii", mfu, U, "U")
display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("u.vtk")
contours = m.contour()
p.add_mesh(m, show_edges=False)
p.add_mesh(contours, color="black", line_width=1)
p.add_mesh(m.contour(8).extract_largest(), opacity=0.1)
pts = m.points
p.show(window_size=[384, 384], cpos="xy")
display.stop()
```

```
[20]:
```

```
<pyvirtualdisplay.display.Display at 0x7f968c03a760>
```

## Exact solution¶

The exact solution to this problem is given by the following equation:

```
[21]:
```

```
evalue = mfu.eval("(1-x*x-y*y)/4")
```

We can calculate the error for the L2 and H1 norms by using `compute`

:

```
[22]:
```

```
L2error = gf.compute(mfu, U - evalue, "L2 norm", mim)
H1error = gf.compute(mfu, U - evalue, "H1 norm", mim)
print("Error in L2 norm : ", L2error)
print("Error in H1 norm : ", H1error)
```

```
Error in L2 norm : 1.965329030567132e-06
Error in H1 norm : 0.00010936971229957788
```

As you can see, the size of the error is within the acceptable range.

## Quadrangular shape case¶

Even if the shape of the region is different, GetFEM can obtain a finite element solution by rewriting the program a little. In the case of a quadrangular shape it may be calculated for:.

```
[23]:
```

```
# Mesh generation
mo = gf.MesherObject("rectangle", [-1.0, -1.0], [1.0, 1.0])
h = 0.1
K = 2
mesh = gf.Mesh("generate", mo, h, K)
# Boundary selection
outer_faces = mesh.outer_faces()
OUTER_BOUND = 1
mesh.set_region(OUTER_BOUND, outer_faces)
sl = gf.Slice(("none",), mesh, 1)
sl.export_to_vtk("sl.vtk", "ascii")
# Mesh draw
display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("sl.vtk")
p.add_mesh(m, show_edges=True)
pts = m.points
p.show(window_size=[512, 384], cpos="xy")
display.stop()
```

```
[23]:
```

```
<pyvirtualdisplay.display.Display at 0x7f96ba6dd040>
```

```
[24]:
```

```
# Definition of finite element methods and integration method
mfu = gf.MeshFem(mesh, 1)
elements_degree = 2
mfu.set_classical_fem(elements_degree)
mim = gf.MeshIm(mesh, pow(elements_degree, 2))
# Model definition
md = gf.Model("real")
md.add_fem_variable("u", mfu)
# Poisson’s equation
md.add_Laplacian_brick(mim, "u")
F = 1.0
md.add_fem_data("F", mfu)
md.set_variable("F", np.repeat(F, mfu.nbdof()))
md.add_source_term_brick(mim, "u", "F")
md.add_Dirichlet_condition_with_multipliers(mim, "u", elements_degree - 1, OUTER_BOUND)
# Model solve
md.solve()
# Export/visualization of the solution
U = md.variable("u")
sl.export_to_vtk("u.vtk", "ascii", mfu, U, "U")
display = Display(visible=0, size=(1280, 1024))
display.start()
p = pv.Plotter()
m = pv.read("u.vtk")
contours = m.contour()
p.add_mesh(m, show_edges=False)
p.add_mesh(contours, color="black", line_width=1)
p.add_mesh(m.contour(8).extract_largest(), opacity=0.1)
pts = m.points
p.show(window_size=[384, 384], cpos="xy")
display.stop()
```

```
[24]:
```

```
<pyvirtualdisplay.display.Display at 0x7f9694d31730>
```

```
[ ]:
```

```
```