# 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.