RobinBCs.py

Contents

  1. Introduction
  2. Variational Formulation
  3. Commented Program
  4. Results
  5. Complete Source Code

Introduction

In this tutorial the implementation of Robin boundary conditions using trace spaces is shown.

The equation solved is the reaction-diffusion equation

\[ - \Delta u + u = 0 \]

in the unit square $ \Omega = (0,1)^2 $ with inhomogeneous Neumann boundary condition

\[ \nabla u\cdot \underline{n} = g(x,y), \qquad (x,y)\in\Gamma_{\rm N}=\{(x,y)\in\partial\Omega\mid y=0\}, \]

and Robin boundary condition

\[ \nabla u\cdot \underline{n} = h_1(x,y) u + h_2(x,y), \qquad (x,y)\in\Gamma_{\rm R}=\partial\Omega\setminus\Gamma_{\rm N} \]

where $ g(x,y) = \sin \pi x $, $ h_1(x,y) = -y $ and $ h_2(x,y) = -\sin \frac{\pi}{2}y $.

Variational Formulation

Now we derive the corresponding variational formulation of the introduced problem: find $ u\in H^1(\Omega) $ such that

\[ \int_{\Omega} \nabla u \cdot \nabla v \;{\rm d}x + \int_{\Omega} u v \;{\rm d}x - \int_{\Gamma_{\rm R}} h_1 u v \;{\rm d}x = \int_{\Gamma_{\rm N}} g v \;{\rm d}x + \int_{\Gamma_{\rm R}} h_2 v \;{\rm d}x \qquad\forall v\in H^1(\Omega). \]

Commented Program

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

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

In order to use MUMPS as a solver, we need to import MPI

10 from mpi4py import MPI

Main Program

The mesh is read from three files containing the coordinates, the elements and the boundary attributes of the mesh.

16  try:
17  mesh = concepts.Import2dMesh(coord = SOURCEDIR + "/applications/unitsquareCoordinates.dat",
18  elms = SOURCEDIR + "/applications/unitsquareElements.dat",
19  boundary = SOURCEDIR + "/applications/unitsquareEdges.dat"
20  )
21  print ('Mesh:')
22  print (mesh)

In our example, the edges are given the following attributes: 1 (bottom), 2 (right), 3 (top) and 4 (left).

If there is an error reading those files, then an exception error message is returned.

24  except RuntimeError:
25  print("Mesh import error")

else, the mesh is plotted using scaling factor of 100, a greyscale of 1.0 and one point per edge.

28  graphics.MeshEPS_r(msh = mesh, filename = "mesh.eps", scale=100, greyscale=1.0, nPoints=1)

Since, we do not have any homogeneous Dirichlet boundary conditions we do not need to set them using bc = concepts.BoundaryConditions().

The formula of the inhomogeneous Neumann data is defined

32  NeumannFormula = concepts.ParsedFormula_r(formula = "(sin(pi*x))")

and a set of unsigned intergers identifying the edge with attribute 1, on which the Neumann boundary condition is prescribed, is created.

34  attrNeumann = concepts.Set_uint(1)

Now, the same is done for the Robin data. The two formulas are defined.

38  RobinFormula1 = concepts.ParsedFormula_r(formula = "(-y)")
39  RobinFormula2 = concepts.ParsedFormula_r(formula = "(-sin(pi/2*y))")

and a set of unsigned intergers identifying the edge with attributes 2, 3 and 4, on which the Robin boundary condition is prescribed, is created.

41  attrRobin = concepts.Set_uint(2)
42  attrRobin.insert(3)
43  attrRobin.insert(4)

In this example, there are no homogeneous Dirichlet boundary conditions and thus, the space can be built without passing any boundary data. We refine the space two times and set the polynomial degree to three. Then, the elements of the space are built and the space is plotted.

47  levelOfRefinement = 2
48  polynomialDegree = 3
49  space = hp2D.hpAdaptiveSpaceH1(msh = mesh, l = levelOfRefinement,
50  p = polynomialDegree)
51  space.rebuild()
52  graphics.MeshEPS_r(spc = space, filename = "space.eps", scale=100, greyscale=1.0, nPoints=1)

Now, the Neumann and Robin trace spaces are built using the 2D space space and the sets attrNeumann and attrRobin of edges with Neumann boundary condition and Robin boundary condition respectively.

56  tscpNeumann = hp2D.TraceSpace(spc = space, edgeAttr = attrNeumann)
57  tscpRobin = hp2D.TraceSpace(spc = space, edgeAttr = attrRobin)

The right hand side is computed. In our case, the vector rhs only contains the integrals of the Neumann and Robin trace as the source term is zero.

61  lformNeumann = hp1D.Riesz_r(frm = NeumannFormula)
62  rhsNeumann = concepts.Vector_r(spc = tscpNeumann, lf = lformNeumann)
63  lformRobin = hp1D.Riesz_r(frm = RobinFormula2)
64  rhsRobin = concepts.Vector_r(spc = tscpRobin, lf = lformRobin)
65  rhs = concepts.Vector_r()
66  rhs = rhsNeumann + rhsRobin;
67 
68  print('\nRHS Vector:')
69  print(rhs)

The system matrix is computed from the bilinear form hp2D.Laplace_r and the bilinear form hp2D.Identity_r.

74  la = hp2D.Laplace_r()
75  A = concepts.SparseMatrix_r(spc = space, bf = la)
76  A.compress()
77 
78  id = hp2D.Identity_r()
79  M = concepts.SparseMatrix_r(spc = space, bf = id)
80  M.compress()
81  M.addInto(A, 1.0)

Furthermore, the integrals of the Robin trace are added.

84  bilformRobin = hp1D.Identity_r(frm = RobinFormula1)
85  M1D = concepts.SparseMatrix_r(spc = tscpRobin, bf = bilformRobin)
86  M1D.compress()
87  M1D.addInto(A, -1.0)
88 
89  print('\nSystem Matrix:')
90  print(A)

We solve the equation using a Mumps solver.

94  solver = concepts.Mumps_r(A = A)
95  sol = concepts.Vector_r(spc = space)
96  solver(fncY = rhs, fncX = sol)
97  print('\nSolver:')
98  print(solver)
99  print('\nSolution:')
100  print(sol)

Finally, the solution u is stored in 'Robin.mat' file using graphics.MatlabBinaryGraphics.To this end, the shape functions are computed on equidistant points using the trapezoidal quadrature rule.

104  hp2D.IntegrableQuad.setTensor(concepts.TRAPEZE, True, 8 )
105  space.recomputeShapefunctions()
106  graphics.MatlabBinaryGraphics(spc = space, filename = "Robin", sol = sol)

Results

Output of the program:

Mesh: Import2dMesh(ncell = 1)
Mesh:
Import2dMesh(ncell = 1)
RHS Vector:
Vector(169, [ 1.549354e-02, 1.678744e-01, 0.000000e+00, -9.444769e-02, 1.570060e-02, 5.980456e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -8.097468e-03, -3.200783e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 2.374103e-01, 0.000000e+00, 3.790460e-02, 2.477186e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -1.745166e-01, 0.000000e+00, 0.000000e+00, -2.305964e-02, -2.713493e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.678744e-01, 0.000000e+00, 3.790460e-02, -2.477186e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 1.549354e-02, -9.444769e-02, 1.570060e-02, -5.980456e-03, -8.097468e-03, -3.200783e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -1.745166e-01, 0.000000e+00, -2.305964e-02, -2.713493e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.280169e-01, -3.451119e-02, -1.813098e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.484019e-01, -2.500000e-01, -4.070872e-02, -6.366754e-04, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.500000e-01, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.280169e-01, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -3.451119e-02, -1.813098e-03, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.500000e-01, -4.166667e-02, 1.734723e-18, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00, -2.484019e-01, -4.166667e-02, 1.734723e-18, -4.070872e-02, -6.366754e-04, 0.000000e+00, 0.000000e+00, 0.000000e+00, 0.000000e+00])
System Matrix:
SparseMatrix(169x169, HashedSparseMatrix: 2825 (9.89111%) entries bound.)
Solver:
Mumps(n = 169)
Solution:
Vector(169, [-2.358621e-01, -1.292237e-01, -2.897389e-01, -3.566355e-01, -4.559817e-02, 1.946218e-02, -1.205068e-02, 1.041054e-03, 7.644850e-03, 1.495471e-03, 6.444423e-02, -1.673132e-02, 1.587023e-03, 1.516002e-02, -1.736858e-02, 2.606585e-04, -6.104195e-02, -2.594242e-01, 6.267945e-02, 3.415858e-03, -4.265729e-02, 3.939918e-03, 2.834458e-02, -1.218636e-03, -3.113952e-02, 3.437624e-03, 4.741337e-04, -5.281713e-04, -4.024244e-01, -4.285123e-01, -1.755653e-02, 8.674954e-04, 2.607856e-02, -7.279306e-06, -1.171772e-02, 3.217059e-04, -4.662949e-03, 1.143960e-03, -7.680474e-04, 4.177298e-05, -5.093568e-01, 3.010143e-02, -9.179211e-04, -1.329014e-02, -1.956168e-03, 1.487323e-02, -1.463781e-03, -2.474342e-03, 6.503121e-04, -1.292237e-01, -2.897389e-01, 6.267945e-02, -3.415858e-03, -1.205068e-02, 1.041054e-03, 2.834458e-02, -1.218636e-03, -3.113952e-02, 3.437624e-03, -4.741337e-04, 5.281713e-04, -2.358621e-01, -3.566355e-01, -4.559817e-02, -1.946218e-02, 6.444423e-02, -1.673132e-02, 7.644850e-03, 1.495471e-03, 1.587023e-03, 1.516002e-02, 1.736858e-02, -2.606585e-04, -5.093568e-01, -4.285123e-01, -1.329014e-02, -1.956168e-03, 3.010143e-02, -9.179211e-04, -1.171772e-02, 3.217059e-04, 1.487323e-02, -1.463781e-03, 2.474342e-03, -6.503121e-04, 2.607856e-02, -7.279306e-06, -4.662949e-03, 1.143960e-03, 7.680474e-04, -4.177298e-05, -5.422839e-01, -5.175787e-01, -1.233352e-02, 2.661886e-04, 2.502922e-02, 2.062011e-04, -1.092909e-02, -4.398691e-04, 1.899724e-03, 1.311650e-04, 3.106070e-04, -2.525089e-05, -6.223890e-01, -2.309497e-02, 9.162962e-05, 3.166459e-02, -1.025598e-03, 7.912876e-03, -6.024182e-05, 5.596927e-04, 5.688176e-05, -6.986325e-01, -6.372163e-01, -1.000702e-02, 2.815548e-03, 1.919949e-02, 1.072334e-03, -5.662442e-03, -7.820322e-04, 4.555847e-03, -2.151253e-03, 2.307172e-03, 1.385236e-05, -6.168261e-01, 2.056321e-02, -1.149454e-04, -4.699151e-03, -6.026103e-04, 8.769205e-04, -2.662965e-04, -8.598998e-05, -8.702793e-05, -5.422839e-01, -6.223890e-01, -1.233352e-02, 2.661886e-04, 3.166459e-02, -1.025598e-03, -2.309497e-02, 9.162962e-05, 7.912876e-03, -6.024182e-05, -5.596927e-04, -5.688176e-05, 2.502922e-02, 2.062011e-04, 1.899724e-03, 1.311650e-04, -3.106070e-04, 2.525089e-05, -6.372163e-01, 2.056321e-02, 1.149454e-04, -5.662442e-03, -7.820322e-04, 8.769205e-04, -2.662965e-04, 8.598998e-05, 8.702793e-05, -6.986325e-01, 1.919949e-02, -1.072334e-03, -1.000702e-02, 2.815548e-03, 4.555847e-03, -2.151253e-03, -2.307172e-03, -1.385236e-05])

Matlab plot of the solution: Matlab output

@section complete Complete Source Code
1 # import modules from package
2 from libconceptspy import concepts
3 from libconceptspy import graphics
4 from libconceptspy import hp1D
5 from libconceptspy import hp2D
6 
7 import os
8 from sys import exit
9 from mpi4py import MPI
10 
11 SOURCEDIR = os.environ['CONCEPTSPATH']
12 
13 def main():
14  try:
15  mesh = concepts.Import2dMesh(coord = SOURCEDIR + "/applications/unitsquareCoordinates.dat",
16  elms = SOURCEDIR + "/applications/unitsquareElements.dat",
17  boundary = SOURCEDIR + "/applications/unitsquareEdges.dat"
18  )
19  print ('Mesh:')
20  print (mesh)
21  except RuntimeError:
22  print("Mesh import error")
23  else:
24  graphics.MeshEPS_r(msh = mesh, filename = "mesh.eps", scale=100, greyscale=1.0, nPoints=1)
25 
26 
27  NeumannFormula = concepts.ParsedFormula_r(formula = "(sin(pi*x))")
28  attrNeumann = concepts.Set_uint(1)
29 
30 
31  RobinFormula1 = concepts.ParsedFormula_r(formula = "(-y)")
32  RobinFormula2 = concepts.ParsedFormula_r(formula = "(-sin(pi/2*y))")
33  attrRobin = concepts.Set_uint(2)
34  attrRobin.insert(3)
35  attrRobin.insert(4)
36 
37 
38  levelOfRefinement = 2
39  polynomialDegree = 3
40  space = hp2D.hpAdaptiveSpaceH1(msh = mesh, l = levelOfRefinement,
41  p = polynomialDegree)
42  space.rebuild()
43  graphics.MeshEPS_r(spc = space, filename = "space.eps", scale=100, greyscale=1.0, nPoints=1)
44 
45 
46  tscpNeumann = hp2D.TraceSpace(spc = space, edgeAttr = attrNeumann)
47  tscpRobin = hp2D.TraceSpace(spc = space, edgeAttr = attrRobin)
48 
49 
50  lformNeumann = hp1D.Riesz_r(frm = NeumannFormula)
51  rhsNeumann = concepts.Vector_r(spc = tscpNeumann, lf = lformNeumann)
52  lformRobin = hp1D.Riesz_r(frm = RobinFormula2)
53  rhsRobin = concepts.Vector_r(spc = tscpRobin, lf = lformRobin)
54  rhs = concepts.Vector_r()
55  rhs = rhsNeumann + rhsRobin;
56 
57  print('\nRHS Vector:')
58  print(rhs)
59 
60 
61 
62  la = hp2D.Laplace_r()
63  A = concepts.SparseMatrix_r(spc = space, bf = la)
64  A.compress()
65 
66  id = hp2D.Identity_r()
67  M = concepts.SparseMatrix_r(spc = space, bf = id)
68  M.compress()
69  M.addInto(A, 1.0)
70 
71  bilformRobin = hp1D.Identity_r(frm = RobinFormula1)
72  M1D = concepts.SparseMatrix_r(spc = tscpRobin, bf = bilformRobin)
73  M1D.compress()
74  M1D.addInto(A, -1.0)
75 
76  print('\nSystem Matrix:')
77  print(A)
78 
79 
80  solver = concepts.Mumps_r(A = A)
81  sol = concepts.Vector_r(spc = space)
82  solver(fncY = rhs, fncX = sol)
83  print('\nSolver:')
84  print(solver)
85  print('\nSolution:')
86  print(sol)
87 
88 
89  hp2D.IntegrableQuad.setTensor(concepts.TRAPEZE, True, 8 )
90  space.recomputeShapefunctions()
91  graphics.MatlabBinaryGraphics(spc = space, filename = "Robin", sol = sol)
92 
93 if __name__ == '__main__':
94  main()
hp2D::TraceSpace
Builds the trace space of an FE space.
Definition: traces.hh:52
graphics::MatlabBinaryGraphics
Class that allows to store graphical infomations in .mat files to use them in Matlab.
Definition: matlabBinaryGraphics.hh:28
Page URL: http://wiki.math.ethz.ch/bin/view/Concepts/WebHome
21 August 2020
© 2020 Eidgenössische Technische Hochschule Zürich