Home » Software » FEATFLOW2 » Concepts » Boundary conditions

Boundary conditions

The topic of boundary conditions is one of the most complicated ones when dealing with partial differential equations. Usually, the choice of the boundary conditions depend on the type of the equation, and thus, the implementation is equation dependent. However, the implementation of boundary conditions is often very similar, and FEAT2 special technique for some of the most common boundary conditions.

Dirichlet boundary conditions

Dirichlet boundary conditions are the most commonly used type of boundary conditions. This type enforces a special value of a function on the boundary. One may think of, e.g., a condition like $$u_{|\partial\Omega}=0$$ for the boundary $\partial\Omega$ of a domain $\Omega$.

There are basically two basic ideas how to implement such boundary conditions:

Overwrite degrees of freedom with the desired value ("strong" implementation)
Introduce penalty terms as boundary integrals ("weak" implementation)

FEAT2 supports the "strong" implementation method (method 1) out-of-the-box. Dirichlet boundary conditions are always "discretised" into a boundary condition structure t_discreteBC (from the module discretebc.f90) which realises nothing else than collecting all degrees of freedom on the boundary together with their corresponding values. In a second step, the t_discreteBC structure is applied to matrices and vectors to modify their content.

A fundamental example

Let us at first give a very basic example to demonstrate the logic behind boundary conditions. We start with the 1D Poisson equation $$ -\partial_{xx} u = f \qquad\text{in $[0,1]$},$$ $$u(0)=u(1)=0$$ for an arbitrary right-hand side $f$. The equation is discretised with the $Q_1$ element on a mesh with two cells and three mesh points: $x_0=0$, $x_1=1/3$, $x_2=2/3$, $x_3=1$, see the following picture.

A discretisation with a linear finite element leads to the following (indefinite) linear system for the unknowns $u_0=u_h(x_0)$, $u_1=u_h(x_1)$, $u_2=u_h(x_2)$, $u_3=u_h(x_3)$:

$$ \begin{pmatrix} 2 & -1 & & \\ -1 & 2 & -1 & \\ & -1 & 2 & -1 \\ & & -1 & 2 \\ \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \end{pmatrix} = \begin{pmatrix} f_0 \\ f_1 \\ f_2 \\ f_3 \\ \end{pmatrix} $$

In the following, we want to implement the Dirichlet boundary conditions

$$ u_0 = 0,\ u_3 = 0.$$

The "strong" implementation of these conditions is done in two steps:

Replace the rows 1 and 4 (corresponding to $x_0$ and $x_3$) by unit vectors.
Replace the right-hand side values $f_0$ and $f_3$ by $u(x_0)=0$ and $u(x_3)=0$.

One obtains the system

$$ \begin{pmatrix} 1 & & & \\ -1 & 2 & -1 & \\ & -1 & 2 & -1 \\ & & & 1 \\ \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \end{pmatrix} = \begin{pmatrix} 0 \\ f_1 \\ f_2 \\ 0 \end{pmatrix}. $$

Solving this system leads to a finite element function $u_h$ which exactly fulfils the desired boundary conditions $u_h(x_0)=u_0=u(x_0)=0$ and $u_h(x_3)=u_3=u(x_3)=0$.

An iterative solution algorithm then sets up a defect in the form $d:=b-Ax$, i.e.,

$$ \begin{pmatrix} d_0 \\ d_1 \\ d_2 \\ d_3 \end{pmatrix} := \begin{pmatrix} 0 \\ f_1 \\ f_2 \\ 0 \end{pmatrix} - \begin{pmatrix} 1 & & & \\ -1 & 2 & -1 & \\ & -1 & 2 & -1 \\ & & & 1 \\ \end{pmatrix} \begin{pmatrix} u_0 \\ u_1 \\ u_2 \\ u_3 \end{pmatrix}. $$

Assuming that $u_0=f_0=u(x_0)$ and $u_3=f_3=u(x_3)$ are set in advance, there must necessarily hold $d_0=d_3=0$ for the defect vector $d=(d_0,d_1,d_2,d_3)$.

Strong implementation of Dirichlet boundary conditions

For the strong implementation of Dirichlet boundary conditions, the following modules are involved:

Module	Description
discretebc.f90	Structure encapsuling boundary conditions
bcassemblybase.f90	Basic structures for the assembly of boundary conditions
bcassembly.f90	Assembly routines for boundary conditions
vectorfilters.f90	Implementation of assmebled boundary conditions into vectors
matrixfilters.f90	Implementation of assmebled boundary conditions into matrices

boundary.f90	Description of regions where boundary conditions are to be applied
meshregion.f90	Description of regions where boundary conditions are to be applied

The core structure for the implementation of boundary conditions is t_discreteBC from discretebc.f90. Boundary conditions are always assembled, which means that the degrees of freedom in combination with the values to prescribe are collected and saved for being imposed to matrices vectors.

Example: In the above example, the structure t_discreteBC contains the following description:

For RHS vectors: Set $f_0:=u(x_0)=0$ and $f_3:=u(x_3)=0$.

For matrices: Replace row 1 and 4 by unit vectors.

For solution vectors: Set $u_0:=u(x_0)=0$ and $u_3:=u(x_3)=0$.

For defect vectors $d$: Set $d_0:=0$ and $d_3:=0$.

There are four steps necessary to work with Dirichlet boundary conditions:

Creation of a structure of type t_discreteBC.
Adding all regions where Dirichlet boundary conditions occur. This necessitates a callback function that provides the values on the boundary. For the description regions on the boundary, mesh regions or boundary regions can be utilised.
Imposing the boundary conditions to the RHS vector, to the system matrix/matrices, probably to the solution vector and all defect vectors.
Releasing the structure.

The following commands are relevant:

Subroutine	Description
bcasm_initDiscreteBC, bcasm_releaseDiscreteBC	Initialises / Releases a boundary condition structure.
bcasm_newDirichletBConMR	Assemble Dirichlet boundary conditions. The desired part of the boundary is described by a mesh region `t_meshRegion`.
bcasm_newDirichletBConRealBd	Assemble Dirichlet boundary conditions. The desired part of the boundary is described by a boundary region `t_boundaryRegion`.
vecfil_discreteBCrhs	Impose the boundary conditions into a RHS vector.
vecfil_discreteBCsol	Impose the boundary conditions into a solution vector.
vecfil_discreteBCdef	Impose the boundary conditions into a defect vector.
matfil_discreteBC	Impose the boundary conditions into a matrix.

The assembly routines bcasm_newDirichletBConMR and bcasm_newDirichletBConRealBd are based on a block discretisation structure t_blockDiscretisation that describes the underlying tensor product space of the equation. Furthermore, these routines accept a parameter iequation which identifies the corresponding subspace where to assemble the boundary conditions for. For example, if a discretisation structure rdiscr describes a finite element space $X_h=V_h\times W_h$, the call

  call bcasm_newDirichletBConMR ( &
      rdiscr, 2, rdiscreteBC, rmeshRegion, fgetBdValues)

assembles boundary conditions for the second subspace, i.e., $W_h$ - corresponding to the 2nd subvector in the solution/RHS vectors and 2nd block row in the corresponding matrices.

Example: Strong implementation of Dirichlet-0 boundary conditions on the complete boundary using mesh regions. The following example demonstrates the above approach on a simple example:

At first, we need a callback function that returns values on the boundary. Let us assume, we want to prescribe $u(x)=0$ for all degrees of freedom on the boundary. We create a proper callback function fgetBdValues.
Second, we initialise a boundary condition structure and add the complete boundary as Dirichlet boundary to it via a mesh region.
Third, we apply the discrete boundary conditions to the RHS, the matrix, the initial solution vector (speeds up the computation) and call the linear solver for solvong the system.
Fourth, we release memory.

The code reads as follows:

subroutine fgetBdValues (...,Dcoords,Dvalues,...)
...
  Dvalues(1) = 0.0_DP
end subroutine

  ...
  type(t_blockDiscretisation) :: rdiscr
  type(t_vectorBlock) :: rx, rrhs
  type(t_matrixBlock) :: rmatrix
  type(t_discreteBC) :: rdiscreteBC
  type(t_meshRegion) :: rmeshRegion
  ...

  ! Initialise the BC structure
  call bcasm_initDiscreteBC(rdiscreteBC)

  ! Get a mesh region descibing the complete boundary
  call mshreg_createFromNodalProp( &
      rmeshRegion, rtriangulation,MSHREG_IDX_ALL)

  ! Assemble boundary conditions into rdiscreteBC.
  call bcasm_newDirichletBConMR ( &
      rdiscr, 1, rdiscreteBC, rmeshRegion, fgetBdValues)

  ! Destroy the mesh region
  call mshreg_done(rmeshRegion)

  ! Apply to the RHS and matrix
  call vecfil_discreteBCrhs (rrhs,rdiscreteBC)
  call matfil_discreteBC (rmatrix,rdiscreteBC)

  ! Apply to the solution vector to have a proper initial solution
  call vecfil_discreteBCsol (rx,rdiscreteBC)

  ! Solve the system
  ...

  ! Release boundary conditions
  call bcasm_releaseDiscreteBC (rdiscreteBC)

One recognizes the following parts of the code:

The subroutine
```
subroutine fgetBdValues (...,Dcoords,Dvalues,...)
...
  Dvalues(1) = 0.0_DP
end subroutine
```
is the callback function which returns Dirichlet boundary values. Dcoords contains the coordinates of the points where values are to be returned. The routine has to set Dvalues(1) to the value in the desired point on the boundary.

The code

call bcasm_initDiscreteBC(rdiscreteBC)
...
call bcasm_releaseDiscreteBC (rdiscreteBC)

creates an empty structure rdiscreteBC and releases it after use.

In the call

call mshreg_createFromNodalProp( &
    rmeshRegion, rtriangulation,MSHREG_IDX_ALL)
...
call mshreg_done(rmeshRegion)

a mesh region is created (and destroyed) which contains all vertices/edges/faces on the complete boundary.

The actual assembly of the boundary conditions is done in
```
call bcasm_newDirichletBConMR ( &
    rdiscr, 1, rdiscreteBC, rmeshRegion, fgetBdValues)
```
This subroutine assembles the boundary conditions based on the mesh region rmeshRegion and calls the callback function fgetBdValues for returning the values in the points on the boundary. Note that the "1" in the above command refers to the first subspace of the finite element tensor space described by rdiscr.
After the boundary conditions are assmbled, the code
```
call vecfil_discreteBCrhs (rrhs,rdiscreteBC)
call matfil_discreteBC (rmatrix,rdiscreteBC)
```
imposes them into the RHS vector and into the matrix. The RHS vector is overwritten by the values of the solution in the points on the boundary, while the corresponding rows in the matrix are replaced by unit vectors.
Finally, the underlying linear system is solved. This is done by
```
call vecfil_discreteBCsol (rx,rdiscreteBC)

! Solve the system
...
```
In a first step, vecfil_discreteBCsol is used to impose the Dirichlet boundary values into the solution vector. This step is optional (as the solution will anyway converge to it), but advised as the linear solver usually converges faster. Then, the actual system is solved.

Example: Inhomogeneous Dirichlet boundary conditions. For inhomogeneous boundary conditions, the callback routine fgetBdValues has to be adjusted properly. The following callback function for example imposes the Dirichlet boundary condition $u(x)=x_1 x_2$ for $x\in\partial\Omega$:

subroutine fgetBdValues (...,Dcoords,Dvalues,...)
...
  dx = Dcoords(1)
  dy = Dcoords(2)
  Dvalues(1) = dx*dy
end subroutine

Example: Dirichlet boundary conditions based on boundary regions. For equations in 2D which involve a description of the boundary via a t_boundary structure, there is an alternative way for describing regions on the boundary where Dirichlet boundary conditions should be imposed. THe user can set up an arbitrary boundary region t_boundaryRegion and apply the routine bcasm_newDirichletBConRealBd. The following code demonstrates this approach:

subroutine fgetBdValues (...,dwhere,Dvalues,...)
...
  real(DP) :: dx,dy

  ! Get x/y-coordinates of the point
  call boundary_getCoords(rdiscretisation%p_rboundary, &
      rboundaryRegion%iboundCompIdx, dwhere, dx, dy)

  Dvalues(1) = ...   ! (use dwhere, dx, dy,...)

end subroutine

  ...
  type(t_discreteBC) :: rdiscreteBC
  type(t_boundaryRegion) :: rregion
  ...

  ! Initialise the BC structure
  call bcasm_initDiscreteBC(rdiscreteBC)

  ! Boundary conditions for 1st bd. component, 1st bd. segment.
  call boundary_createRegion (rboundary, 1, 1, rregion)
  call bcasm_newDirichletBConRealBd ( &
      rdiscr, 1, rregion, rdiscreteBC, fgetBdValues)

  ! Boundary conditions for 1st bd. component, 2nd bd. segment.
  call boundary_createRegion (rboundary, 1, 2, rregion)
  call bcasm_newDirichletBConRealBd ( &
      rdiscr, 1, rregion, rdiscreteBC, fgetBdValues)
  ...

In this example, boundary_createRegion from boundary.f90 is invoked to create a boundary region for segment 1, 2, ... on the 1st boundary component of the the analytical boundary. In a second step, bcasm_newDirichletBConRealBd is called to assemble the boundary conditions on each boundary region.

Remark: In contrast to mesh regions, the callback routine fgetBdValues can not directly access the x/y coordinates of the point where the boundary condition is to be imposed. Instead, the point is identified via a parameter value dwhere and a boundary component rboundaryRegion%iboundCompIdx. The parameter value dwhere is given in 0-1 parametrisation. If necessary, one can obtain the corresponding x/y coordinates via the routine boundary_getCoords as indicated in the above example.

Example: Parametrised boundary conditions. Similar to all callback-based routines, the callback routine fgetBdValues can be parametrised via a t_collection structure from collection.f90. For this purpose, the usual "collection strategy" of the FEAT2 package is used. The assembly routine accepts an optional parameter rcollection of type t_collection which is passed to the callback function:

Type of data	passed via
Double precision data can be passed via	`(t_collection)%DquickAccess(:)`
Integers can be passed via	`(t_collection)%IquickAccess(:)`
Strings can be passed via	`(t_collection)%SquickAccess(:)`
Block vectors can be passed via	`(t_collection)%p_rvectorQuickAccessX` with N=1,2,...
Parsers can be passed via	`(t_collection)%p_rfparserQuickAccessN` with N=1,2,...
Parameter lists can be passed via	`(t_collection)%p_rparlistQuickAccessN` with N=1,2,...
Complex parameters can be specified via	named items in the collection.

etc. Using only the quick access arrays, the collection does not even have to be initialised.

The following example demonstrates how to set up a parabolic heating on the 1st boundary segment, where themaximum temperature is specified in a variable dmaxTemp:

subroutine fgetBdValues (...,dwhere,Dvalues,rcollection)
...
  real(DP) :: dmaxTemp
  integer :: iequation

  ! Get the parameters
  dmaxTemp = rcollection%DquickAccess(1)
  iequation = rcollection%IquickAccess(1)

  if (iequation .eq. 1) then
    Dvalues(1) = mprim_getParabolicProfile (dwhere,1.0_DP,dmaxTemp)
  ...
end subroutine

  ...
  type(t_discreteBC) :: rdiscreteBC
  type(t_boundaryRegion) :: rregion

  real(DP) :: dmaxTemp
  type(t_collection) :: rcollection
  ...

  ! Initialise the BC structure
  call bcasm_initDiscreteBC(rdiscreteBC)

  ! Prepare a collection structure, pass dmaxTemp
  ! and the current equation
  rcollection%IquickAccess(1) = 1
  rcollection%DquickAccess(1) = dmaxTemp

  ! Boundary conditions for 1st bd. component, 1st bd. segment.
  call boundary_createRegion (rboundary, 1, 1, rregion)
  call bcasm_newDirichletBConRealBd (rdiscr, rcollection%IquickAccess(1), &
      rregion, rdiscreteBC, fgetBdValues, rcollection)
  ...

The parameters dmaxTemp as well as the number of the current equation are saved to DquickAccess(1) and IquickAccess(1) and passed via the collection to the callback routine. The callback routine can then extract the parameter from the collection. Passing the number of the current equation has been proven useful, so the callback routine can decide what is to be assembled.

Do-nothing (Neumann) boundary conditions

This type of boundary condition is the most simple one. One just does - nothing. A matrix and a right-hand side vector coming from the matrix/vector assembly already "prepared" for Neumann boundary conditions and do not have to be modified by any means to implement this. However, one has to note that pure do-nothing boundary conditions usually lead to an indefinite system, and thus, further modifications often have to be applied to the system to obtain a definite one.

Featflow

Boundary conditions

Dirichlet boundary conditions

A fundamental example

Strong implementation of Dirichlet boundary conditions

Do-nothing (Neumann) boundary conditions