howToGetStarted.py

On this page you find an easy example of an elliptic partial differential equation (PDE) solved with Concepts using Python bindings. Starting with this PDE, we derive its associated variational formulation and solve it numerically. Using the example of the standard concepts setting howToGetStarted.cc in the Concepts tutorials we show the usage of Concepts Python bindings.

Introduction

We want to solve the following elliptic PDE with homogeneous Dirichlet boundary conditions

\begin{eqnarray*} - \Delta u + u &=& f \quad \textrm{in} \quad \Omega\\ u &=& 0 \quad \textrm{on} \quad \partial \Omega\ \end{eqnarray*}

with the right hand side $f(x,y) = (2 \pi^2 +1)\sin(\pi x) \sin(\pi y) $ and $ \Omega = (0,1)^2$. This problem has the unique solution $u(x,y) = \sin(\pi x) \sin(\pi y) $.

To solve this problem with Concepts, we first need its variational formulation. Let $ H_0^1(\Omega)$ be the Sobolev space of square integrable and weak differentiable functions with vanishing trace on $ \partial \Omega$. Then we seek the weak solution $ u \in H_0^1(\Omega) $ that satisfies

\[ - \int_{\Omega} \Delta u v \, {\rm d}\mathbf{x} + \int_\Omega u v \, {\rm d}\mathbf{x} = \int_\Omega f v \, {\rm d}\mathbf{x} \qquad \forall v \in H_0^1(\Omega). \]

Using integration by parts (Gauss theorem), this is equivalent to

\[ \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d}\mathbf{x} - \int_{\partial \Omega} \underline{n} \cdot \nabla u \, v \, {\rm d}s(\mathbf{x})+ \int_\Omega u v \, {\rm d}\mathbf{x} = \int_\Omega f v \, {\rm d}\mathbf{x} \qquad \forall v \in H_0^1(\Omega), \]

where $ \underline{n} $ denotes the unit normal vector outward to $ \Omega $. Using that $ v \in H_0^1(\Omega)$ has vanishing trace, we get

\[ \int_{\Omega} \nabla u \cdot \nabla v \, {\rm d}\mathbf{x} + \int_\Omega u v \, {\rm d}\mathbf{x} = \int_\Omega f v \, {\rm d}\mathbf{x} \qquad \forall v \in H_0^1(\Omega). \]

This is the variational formulation of our problem. Let now $ V_h \subseteq H_0^1$ be a finite dimensional subspace of $ H_0^1$. The Galerkin approximation of our solution $u$ in the space $ V_h$ is given as the (unique) solution $u_h \in V_h$ that fulfills the equation

\[ \int_{\Omega} \nabla u_h \cdot \nabla v_h \, {\rm d}\mathbf{x} + \int_\Omega u_h v_h \, {\rm d}\mathbf{x} = \int_\Omega f v_h {\rm d}\mathbf{x} \qquad \forall v_h \in V_h. \]

Let $ \phi_1, ..., \phi_N$ be a basis of the space $ V_h $, then we can find a unique vector $\alpha = (\alpha_i)_{i = 1,...,N}$ such that $ \sum_i \alpha_i \phi_i = u_h$. Plugging this into the equation above and using the linearity of the integral, we obtain the linear equation

\[ \sum_k \alpha_k \left( \int_\Omega \nabla \phi_k \cdot \nabla \phi_j \, {\rm d}\mathbf{x} + \int_\Omega \phi_k \phi_j \, {\rm d}\mathbf{x} \right) = \int_\Omega f \phi_j \, {\rm d}\mathbf{x} \qquad \forall j \in \{ 1, ..., N\}. \]

For $s_{k,j}:= \int_\Omega \nabla \phi_k \nabla \phi_j \,{\rm d}\mathbf{x} $, $m_{k,j}:= \int_\Omega \phi_k \phi_j \,{\rm d}\mathbf{x} $ and $b_{j}:= \int_\Omega f \phi_j \,{\rm d}\mathbf{x} $ and matrices $ S = (s_{ij})_{i j}$, $ M = (m_{ij})_{i j}$ and a vector $b = (b_i)_i$ this is equivalent to the linear equation

\[ (S + M) \alpha = b \]

with the stiffness matrix $ S $, the mass matrix $ M $ and the load vector $ b $ of the right hand side.

Concepts uses the hp-finite element method to approximate the solution of a PDE, which means it calculates the Galerkin approximation of the solution in the finite dimensional subspace

\[ V_{M_h, p} := \left\{ u \in H^1_0 \textrm{ such that } u\mid_\tau \textrm{ is a polynomial of degree } \leq p \quad \forall \tau \in M_h \right\} \]

where $ M_h $ is a mesh on $ \Omega $ with mesh size $ h $ and $ \tau \in M_h $ is a cell in the mesh.

Commented Program

Now let us see how to calculate the FEM-approximation of the solution of the PDE with Concepts.

First, we have to import certain modules from the libconceptspy package.

2 from libconceptspy import concepts
3 from libconceptspy import hp2D
4 from libconceptspy import graphics

The module concepts contains essential functionalities and must be loaded at first. The module hp2D includes all functions that are needed to build a representation of the hp-finite element space in two dimensions. The module graphics is used for the graphical output of the mesh and the solution.

Next we start with the actual code, consisting of a main function.

10 def main():

First we build the geometry, a square with one vertex in the origin, and another one in the point $ (1,1) $. The third parameter is the attribute of all edges, which we set to to $ 1 $.

14  mesh = concepts.Square(sizex=1, sizey=1, attrib=1)

Now we prescribe homogeneous Dirichlet boundary conditions to all edges of attribute $ 1 $, i.e. to all edges.

17  boundary = concepts.Boundary( type=concepts.Boundary.DIRICHLET )
18  boundaryConditions = concepts.BoundaryConditions()
19  boundaryConditions.add( attrib = 1, bcObject=boundary )

With the mesh and the boundary conditions at hand we can build the hp-finite element space.

22  levelOfRefinement = 2
23  polynomialDegree = 6
24 
25  space = hp2D.hpAdaptiveSpaceH1( msh = mesh, l = levelOfRefinement,
26  p = polynomialDegree, bc = boundaryConditions )
27  space.rebuild()

This is a finite element space on the square, with a twice refined mesh (a single refinement of a two dimensional mesh divides an element in the mesh into four elements, so in our case we have 16 subelements in our square) and polynomial degree six on every element.

Now we define the linear form on the right hand side of the variational formulation.

31  formula = concepts.ParsedFormula_r( formula = "((2*pi*pi+1)*sin(pi*x)*sin(pi*y))" )
32  lform = hp2D.Riesz_r( frm = formula )

In the first line we define a representation of $ f $ and in the second line we define a representation of the linear form $ v \mapsto \int_\Omega f v \, {\rm d}\mathbf{x} $.

In the next line we calculate the vector of the right hand side using the space space and the linear form lform.

35  rhs = concepts.Vector_r( spc = space, lf = lform )

We compute the stiffness and mass matrices in a similar way. We first initialize the bilinear forms.

38  la = hp2D.Laplace_r()
39  id = hp2D.Identity_r()

The first one represents the bilinear form $ (u,v) \mapsto \int_\Omega \nabla u \cdot \nabla v \, {\rm d}\mathbf{x} $ that is the variational representation of the Laplace operator and the second one represents the mapping $ (u,v) \mapsto \int_\Omega u v \, {\rm d}\mathbf{x} $.

Now we build the matrices using a concepts.SparseMatrix_r class that computes the matrix belonging the space space and the bilinear forms la and id, respectively.

42  S = concepts.SparseMatrix_r( spc = space, bf = la)
43  M = concepts.SparseMatrix_r( spc = space, bf = id)

Before we solve the linear equation, we add the mass matrix into the stiffness matrix.

46  M.addInto(S, 1.0)

We solve the linear equation using a conjugate gradient (CG) solver which is appropriate as the the system matrix $ S+M $ is symmetric. We set the error tolerance to 1e-12 and the maximum number of iterations to 400.

53  solver = concepts.CG_r(A = S, maxeps = 1e-6, maxit = 200)
54 
55  sol = concepts.Vector_r(spc = space)
56 
57  solver(fncY = rhs, fncX = sol)

Now, we want to save the information about the mesh and the values of the solution on the integration points in a MAT-file that can be read by Matlab or gnu Octave.

First, we set the integration rule to the trapezoidal rule and recompute the shape functions on this points.

61  hp2D.IntegrableQuad.setTensor( concepts.TRAPEZE, True, polynomialDegree )
62  space.recomputeShapefunctions()

Then we save the information about the mesh in a MAT-file named "matlabSolution.mat" and add the values of the solution on the integration points as a vector named "u".

64  mbg = graphics.MatlabBinaryGraphics(spc = space, filename = "matlabSolution")
65  mbg.addSolution(spc = space, name = "u", sol = sol)

Finally we save the information also in a .vtk-file.

69  vtk_g = graphics.VtkGraphics(spc= space, filename = "vtkSolution")
70  vtk_g.addSolution(spc = space, name = "u", sol = sol)

By running the program you compute a FE-approximation of the solution of the PDE. In the directory where you run the program you will find a MAT-file named "matlabSolution.mat" and a .vtk file named "vtkSolution.vtk". If you use Matlab you can load MAT-file and plot the solution entering the following instructions in your Matlab command line

> load firstSolution.mat
> patch(x(msh), y(msh), u(msh))

This gives you the following visualization of the solution $ u $. Matlab output

Complete Source Code

1 # import modules from package
2 from libconceptspy import concepts
3 from libconceptspy import hp2D
4 from libconceptspy import graphics
5 
6 if concepts.HAS_MPI :
7  from mpi4py import MPI
8 
9 def main():
10 
11 
12  mesh = concepts.Square(sizex=1, sizey=1, attrib=1)
13 
14  boundary = concepts.Boundary( type=concepts.Boundary.DIRICHLET )
15  boundaryConditions = concepts.BoundaryConditions()
16  boundaryConditions.add( attrib = 1, bcObject=boundary )
17 
18  levelOfRefinement = 2
19  polynomialDegree = 6
20 
21  space = hp2D.hpAdaptiveSpaceH1( msh = mesh, l = levelOfRefinement,
22  p = polynomialDegree, bc = boundaryConditions )
23  space.rebuild()
24 
25 
26  formula = concepts.ParsedFormula_r( formula = "((2*pi*pi+1)*sin(pi*x)*sin(pi*y))" )
27  lform = hp2D.Riesz_r( frm = formula )
28 
29  rhs = concepts.Vector_r( spc = space, lf = lform )
30 
31  la = hp2D.Laplace_r()
32  id = hp2D.Identity_r()
33 
34  S = concepts.SparseMatrix_r( spc = space, bf = la)
35  M = concepts.SparseMatrix_r( spc = space, bf = id)
36 
37  M.addInto(S, 1.0)
38 
39  # Solver
40  if concepts.HAS_MUMPS :
41  solver = concepts.Mumps_r(A = S)
42  else :
43  solver = concepts.CG_r(A = S, maxeps = 1e-6, maxit = 200)
44 
45  sol = concepts.Vector_r(spc = space)
46 
47  solver(fncY = rhs, fncX = sol)
48 
49 
50  hp2D.IntegrableQuad.setTensor( concepts.TRAPEZE, True, polynomialDegree )
51  space.recomputeShapefunctions()
52  mbg = graphics.MatlabBinaryGraphics(spc = space, filename = "matlabSolution")
53  mbg.addSolution(spc = space, name = "u", sol = sol)
54 
55 
56  vtk_g = graphics.VtkGraphics(spc= space, filename = "vtkSolution")
57  vtk_g.addSolution(spc = space, name = "u", sol = sol)
58 
59 if __name__== "__main__" :
60  main()
graphics::MatlabBinaryGraphics
Class that allows to store graphical infomations in .mat files to use them in Matlab.
Definition: matlabBinaryGraphics.hh:28
concepts::Boundary
Class to describe an element of the boundary.
Definition: boundary.hh:35
concepts::Square
Mesh for with one quadrilateral.
Definition: square.hh:22
graphics::VtkGraphics
Class that allows to store graphical infomations in .vtk files to use them in paraview.
Definition: vtkGraphics.hh:20
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich