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.
from IPython.display import YouTubeVideo YouTubeVideo("6JuB1GiDLQQ", start=512)
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 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.
GetFEM can be imported as follows (numpy has also to be imported).
import getfem as gf import numpy as np
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.
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).
h = 0.1 K = 2 mesh = gf.Mesh("generate", mo, h, K)
To define a boundary condition, we set a boundary number on the outer circumference of the circle.
outer_faces = mesh.outer_faces() OUTER_BOUND = 1 mesh.set_region(OUTER_BOUND, outer_faces)
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.
sl = gf.Slice(("none",), mesh, 1)
We can export a slice to VTK file by using export_to_vtk method.
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.
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()
<pyvirtualdisplay.display.Display at 0x7f6be0527c40>
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).
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.
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.
mim = gf.MeshIm(mesh, pow(elements_degree, 2))
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).
md = gf.Model("real") md.add_fem_variable("u", mfu)
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
If you want to define constants in GetFEM, we use the
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.
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.
md.add_source_term_brick(mim, "u", "F")
Finally, we set the Dirichlet condition at the boundary.
md.add_Dirichlet_condition_with_multipliers(mim, "u", elements_degree - 1, OUTER_BOUND)
Once the model is correctly defined, we can simply solve it by:
Export/visualization of the solution¶
The finite element problem is now solved. We can get the solution \(u\) by using
U = md.variable("u")
We can output the computed u with the mesh of the Slice object.
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") p.add_mesh(m, show_edges=True) pts = m.points p.show(window_size=[512, 384], cpos="xy") display.stop()
<pyvirtualdisplay.display.Display at 0x7f6b9815f1c0>
The exact solution to this problem is given by the following equation:
evalue = mfu.eval("(1-x*x-y*y)/4")
We can calculate the error for the L2 and H1 norms by using
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.