This tutorial goes in depth into the mechanisms required to solve the Dirichlet problem
-\Delta u=f \ \ \ in \ \ \Omega
with a nonzero Dirichlet boundary condition
u|_{\Gamma_D}=g \ \ \ on \ a \ boundary \ part \ \Gamma_D \subset \partial \Omega
The same mechanisms are used in solving boundary value problems involving operators other than the Laplacian.
You will see how to perform these tasks in NGSolve:
- extend Dirichlet data from boundary parts,
- convert boundary data into a volume source,
- reduce inhomogeneous Dirichlet case to the homogeneous case, and
- perform all these tasks automatically within a utility.
Spaces with Dirichlet conditions on part of the boundary
from ngsolve import *
from ngsolve.webgui import Draw
from netgen.geom2d import unit_square
mesh = Mesh(unit_square.GenerateMesh(maxh=0.2))
mesh.GetBoundaries()
#Return
('bottom', 'right', 'top', 'left')
The first line of code you provided creates a mesh object from the unit_square
object, and the GenerateMesh
method is used to generate a mesh with a maximum element size (or “maximum h-size”) of 0.2 units. The second line of code mesh.GetBoundaries()
is used to retrieve the boundary of the mesh. It will return the boundary information of the mesh, like boundary ID, boundary name, as well as boundary elements as a tuple of the element index and face index. This can be useful in setting up boundary conditions in a numerical simulation.
Please do note that this is a 2D example and if you want to work with 3D meshes you can use unit_cube or unit_tet object to generate 3D meshes.
The unit_square has its boundaries marked as left, right, top and bottom. Suppose we want non-homogeneous Dirichlet boundary conditions on
\Gamma_D = \Gamma_{left} \cup \Gamma_{right}
Then, we set the space as follows:
fes = H1(mesh, order=2, dirichlet="left|right")
The line of code you provided creates a finite element space (FES) object, which is used to represent the function space in which the solution to a partial differential equation (PDE) will be sought.
The H1
function creates a finite element space that consists of continuous, piecewise polynomials of a specified order. In this case, the order of the polynomials is 2. This means that the elements of the FES will be composed of polynomials of degree 2 or less, defined over the mesh.
The mesh
object that is passed as an argument is the mesh over which the FES is defined.
The dirichlet
argument is used to specify which boundaries of the mesh will have a Dirichlet boundary condition applied. In this case, “left” and “right” are specified, which means that the left and right boundaries of the mesh will have Dirichlet boundary conditions applied.
This FES object can now be used to perform numerical computations using the finite element method, such as solving a PDE or computing the energy of a given solution.
Compare this space with the one without the dirichlet
flag:
fs2 = H1(mesh, order=2)
fes.ndof, fs2.ndof # total number of unknowns
#return
(125, 125)
This code creates another FES object fs2
using the same mesh, but with a different configuration. In this case, no dirichlet
argument is provided, meaning that no Dirichlet boundary condition is imposed on this FES.
The ndof
attribute of the FES objects fes
and fs2
returns the total number of degrees of freedom (DOF) in the finite element space. In other words, it gives the number of unknowns in the system of equations that will be solved to find the solution of a PDE.
fes.ndof and fs2.ndof are the total number of DOF of fes and fs2 FES, respectively. Since fes and fs2 have the same order and defined on the same mesh, the difference between fes.ndof and fs2.ndof is the number of degree of freedom on the left and right bounday which were constrained. If you want to check for the specific number of dof on the boundaries, you can use fes.GetDofs(BoundaryFlags) function.
Thus, the dirichlet
flag did not change ndof
. In NGSolve the unknowns are split into two groups: * dirichlet dofs (or constrained dofs), * free dofs.
The facility FreeDofs
gives a BitArray
such that FreeDofs[dof]
is True if and only if dof
is a free degree of freedom.
print("free dofs of fs2 without \"dirichlet\" flag:\n",
fs2.FreeDofs())
print("free dofs of fes:\n", fes.FreeDofs())
#return
free dofs of fs2 without "dirichlet" flag:
0: 11111111111111111111111111111111111111111111111111
50: 11111111111111111111111111111111111111111111111111
100: 1111111111111111111111111
free dofs of fes:
0: 00001111000011110000111111111111111111010011011111
50: 11111110110110111111111111111010111011111111111111
100: 1111111111111111111111111
The FreeDofs()
method returns a list of indices corresponding to the degrees of freedom that are not constrained by any condition (i.e., not fixed to a particular value). These are the DOF that are solved for when a PDE is solved using the finite element method.
In the first line, it prints the list of free dofs of fs2 FES which don’t have any “dirichlet” flag and in the second line, it prints the free dofs of fes FES.
As you already know that fes object have “dirichlet” flag on left and right side. Therefore, fes.FreeDofs()
will return a list of indices corresponding to the degrees of freedom which are not constrained by any condition, except for the left and right boundaries. While fs2.FreeDofs()
will return a list of indices corresponding to the degrees of freedom which are not constrained by any condition.
By comparing the free dofs of fs2 and fes, you can see how many degree of freedom are fixed by the dirichlet flag in the fes FES.
Leave a Reply