Boundary Conditions
$ \newcommand{\cross}{\times} \newcommand{\inner}{\cdot} \newcommand{\div}{\nabla\cdot} \newcommand{\curl}{\nabla\times} \newcommand{\grad}{\nabla} \newcommand{\ddx}[1]{\frac{d#1}{dx}} \newcommand{\abs}[1]{|#1|} \newcommand{\dO}{{\partial\Omega}} $
MFEM supports boundary conditions of mixed type through the definition of boundary attributes on the mesh. A boundary attribute is a positive integer assigned to each boundary element of the mesh. Since each boundary element can have only one attribute number the boundary attributes split the boundary into a group of disjoint sets. MFEM allows the user to define boundary conditions on a subset of boundary attributes.
Typically mixed boundary conditions are imposed on disjoint portions of the boundary defined as:
Symbol | Description |
---|---|
$\Gamma\equiv\dO$ | Boundary of the Domain ($\Omega$) |
$\Gamma_D$ | Dirichlet Boundary |
$\Gamma_N$ | Neumann Boundary |
$\Gamma_R$ | Robin Boundary |
$\Gamma_0$ | Natural Boundary |
Where we assume $\Gamma = \Gamma_D\cup\Gamma_N\cup\Gamma_R\cup\Gamma_0$. In MFEM boundaries are usually described by "marker arrays". A marker array is an array of integers containing zeros and ones with a length equal to the largest boundary attribute index.
// Assume we start with an array containing boundary attribute numbers
// stored in bdr_attr.
//
// Prepare a marker array from a set of attributes
Array<int> bdr_marker(pmesh.bdr_attributes.Max());
bdr_marker = 0;
for (int i=0; i<bdr_attr.Size(); i++)
{
bdr_marker[bdr_attr[i]-1] = 1;
}
Separate marker arrays of this type can be prepared for the Dirichlet, Neumann, and Robin portions of the boundary. It is left to the user to ensure that $\Gamma_D$, $\Gamma_N$, and $\Gamma_R$ are defined as non-overlapping portions of the boundary.
Continuous Formulations
Dirichlet (Essential) Boundary Conditions
In continuous formulations essential boundary conditions are set by modifying the linear system to require the degrees of freedom on the boundary to obtain specific values. This limits the types of constraints that can be imposed on fields. For example, $L^2$ fields have no degrees of freedom on the boundary of elements so essential BCs cannot be applied, H(Curl) (Nedelec) elements can only constrain the tangential components of a vector field, and H(Div) (Raviart-Thomas) elements can only constrain the normal component of a vector field.
Space | Essential BC |
---|---|
H1 | $u = f$ on $\Gamma_D$ |
H1$^d$ | $\vec{u} = \vec{f}$ on $\Gamma_D$ |
ND | $(\hat{n}\times\vec{u})\times\hat{n} = \vec{f}$ on $\Gamma_D$ |
RT | $\hat{n}\cdot\vec{u} = f$ on $\Gamma_D$ |
MFEM provides a convenience method, called FormLinearSystem
, on the
[Par]BilinearForm
class which can prepare a linear system with these
essential constraints.
// Set the Dirichlet values in the solution vector
ParGridFunction u(&fespace);
u = 0.0;
u.ProjectBdrCoefficient(dbcCoef, dbc_marker);
// Prepare the source term in the right-hand-side
ParLinearForm b(&fespace);
b.AddDomainIntegrator(new DomainLFIntegrator(rhsCoef));
b.Assemble();
// Prepare the bilinear form
ParBilinearForm a(&fespace);
a.AddDomainIntegrator(new DiffusionIntegrator(matCoef));
a.Assemble();
// Determine the essential degrees of freedom corresponding to the set of
// boundary attributes marked in dbc_marker
Array<int> ess_tdof_list(0);
fespace.GetEssentialTrueDofs(dbc_marker, ess_tdof_list);
// Prepare the linear system with enforcement of the essential boundary
// conditions
OperatorPtr A;
Vector B, X;
a.FormLinearSystem(ess_tdof_list, u, b, A, X, B);
Natural Boundary Conditions
The so called "Natural Boundary Conditions" arise whenever weak derivatives occur in a PDE (see below for more on weak derivatives). Weak derivatives must be handled using integration by parts which introduces a boundary integral. If this boundary integral is ignored, its value is implicitly set to zero which creates an implicit constraint on the solution called a "natural boundary condition".
Continuous Operator | Weak Operator | Natural BC |
---|---|---|
$-\div(\lambda\grad u)$ | $(\lambda\grad u,\grad v)$ | $\hat{n}\cdot(\lambda\grad u)=0$ on $\Gamma_0$ |
$\curl(\lambda\curl\vec{u})$ | $(\lambda\curl\vec{u},\curl\vec{v})$ | $\hat{n}\cross(\lambda\curl\vec{u})=0$ on $\Gamma_0$ |
$-\grad(\lambda\div\vec{u})$ | $(\lambda\div\vec{u},\div\vec{v})$ | $\lambda\div\vec{u}=0$ on $\Gamma_0$ |
$\div(\vec{\lambda}u)$ | $(-\vec{\lambda}u,\grad v)$ | $\hat{n}\cdot\vec{\lambda}u = 0$ on $\Gamma_0$ |
$\curl(\lambda\vec{u})$ | $(\lambda\vec{u},\curl\vec{v})$ | $\hat{n}\cross(\lambda\vec{u})=0$ on $\Gamma_0$ |
$-\div(\lambda\grad u) + \div(\vec{\beta}u)$ | $(\lambda\grad u - \vec{\beta}u,\grad v)$ | $\hat{n}\cdot(\lambda\grad u-\vec{\beta}u)=0$ on $\Gamma_0$ |
No additional implementation is necessary to impose natural boundary conditions. Any portion of the boundary where a Dirichlet, Neumann, or Robin boundary condition has not been applied will receive a natural boundary condition by default.
Neumann Boundary Conditions
Neumann boundary conditions are closely related to natural boundary conditions. Rather than ignoring the boundary integral we integrate a known function on the boundary which approximates the desired value of the boundary condition (often a involving a derivative of the field). The following table shows a variety of common operators and their related Neumann boundary condition.
Operator | Continuous Operator | Neumann BC |
---|---|---|
$(\lambda\grad u,\grad v)$ | $-\div(\lambda\grad u)$ | $\hat{n}\cdot(\lambda\grad u)=f$ on $\Gamma_N$ |
$(\lambda\curl\vec{u},\curl\vec{v})$ | $\curl(\lambda\curl\vec{u})$ | $\hat{n}\cross(\lambda\curl\vec{u})=\hat{n}\cross\vec{f}$ on $\Gamma_N$ |
$(\lambda\div\vec{u},\div\vec{v})$ | $-\grad(\lambda\div\vec{u})$ | $\lambda\div\vec{u}=\hat{n}\cdot\vec{f}$ on $\Gamma_N$ |
$(-\vec{\lambda}u,\grad v)$ | $\div(\vec{\lambda}u)$ | $\hat{n}\cdot\vec{\lambda}u = f$ on $\Gamma_N$ |
$(\lambda\vec{u},\curl\vec{v})$ | $\curl(\lambda\vec{u})$ | $\hat{n}\cross(\lambda\vec{u})=\hat{n}\cross\vec{f}$ on $\Gamma_N$ |
$(\lambda\grad u - \vec{\beta}u,\grad v)$ | $-\div(\lambda\grad u) + \div(\vec{\beta}u)$ | $\hat{n}\cdot(\lambda\grad u-\vec{\beta}u)=f$ on $\Gamma_N$ |
To impose these boundary conditions in MFEM simply modify the
right-hand side of your linear system by adding the appropriate
boundary integral of either $f$ or $\vec{f}$. For $H^1$ or $L^2$
fields this can be accomplished by adding the BoundaryLFIntegrator
with an appropriate coefficient for $f$ to a [Par]LinearForm
object.
Neumann boundary conditions can be added to the above example code by adding the following line before the call to b.Assemble()
.
// Add Neumann BCs n.(matCoef Grad u) = nbcCoef on the boundary marked in
// the nbc_marker array.
b.AddBoundaryIntegrator(new BoundaryLFIntegrator(nbcCoef), nbc_marker);
For H(Curl) fields this can be accomplished by adding the
VectorFEBoundaryTangentLFIntegrator
with an appropriate vector
coefficient for $\vec{f}$ to a [Par]LinearForm
object. And
finally, for H(Div) fields this can be accomplished by adding the
VectorFEBoundaryFluxLFIntegrator
with an appropriate scalar
coefficient for $f = \hat{n}\cdot\vec{f}$ to a [Par]LinearForm
object. Other integrators may be appropriate if it is desirable to
express the functions $\,f$ or $\vec{f}$ in other ways.
Robin Boundary Conditions
Robin boundary conditions typically involve a linear function of the field and its normal derivative. As such they also arise from weak derivatives and the boundary integrals they introduce to the system of equations. Typical forms of the Robin boundary condition include the following.
Operator | Continuous Operator | Robin BC |
---|---|---|
$(\lambda\grad u,\grad v)$ | $-\div(\lambda\grad u)$ | $\hat{n}\cdot(\lambda\grad u)+\gamma\,u=f$ on $\Gamma_R$ |
$(\lambda\curl\vec{u},\curl\vec{v})$ | $\curl(\lambda\curl\vec{u})$ | $\hat{n}\cross(\lambda\curl\vec{u}+\gamma\,\hat{n}\cross\vec{u})=\hat{n}\cross\vec{f}$ on $\Gamma_R$ |
$(\lambda\div\vec{u},\div\vec{v})$ | $-\grad(\lambda\div\vec{u})$ | $\lambda\div\vec{u}+\gamma\,\hat{n}\cdot\vec{u}=\hat{n}\cdot\vec{f}$ on $\Gamma_R$ |
$(\lambda\grad u - \vec{\beta}u,\grad v)$ | $-\div(\lambda\grad u) + \div(\vec{\beta}u)$ | $\hat{n}\cdot(\lambda\grad u-\vec{\beta}u)+\gamma\,u=f$ on $\Gamma_R$ |
Robin boundary conditions are applied in the same manner as Neumann
boundary conditions except that one must also add a boundary integral
to the [Par]BilinearForm
object to account for the term involving
$\gamma$. For example, when solving for an $H^1$ field one should add
a MassIntegrator
with an appropriate scalar coefficient for
$\gamma$.
The implementation of a Robin boundary condition requires precisely
the same change to the right-hand-side as the Neumann boundary
condition as well as a new term in the bilinear form before a.Assemble()
:
// Add Robin BCs n.(matCoef Grad u) + rbcACoef u = rbcBCoef on the boundary
// marked in the rbc_marker array.
b.AddBoundaryIntegrator(new BoundaryLFIntegrator(rbcBCoef), rbc_marker);
...
// Add Robin BCs n.(matCoef grad u) + rbcACoef u = rbcBCoef on the boundary
// marked in the rbc_marker array.
a.AddBoundaryIntegrator(new MassIntegrator(rbcACoef), rbc_marker);
Discontinuous Galerkin Formulations
In the Discontinuous Galerkin (DG) formulation the
Natural,
Neumann, and Robin
can be implemented in a similar the same manner as in the continuous case
(adding the appropriate LinearFormIntegrator
as a boundary face integrator
instead of a boundary integrator). However, since DG basis functions
have no degrees of freedom associated with the boundary, Dirichlet boundary
conditions must be handled differently.
// Add the desired value for the Dirichlet constraint on the boundary
// marked in the dbc_marker array.
b.AddBdrFaceIntegrator(new DGDirichletLFIntegrator(dbcCoef, matCoef, sigma, kappa),
dbc_marker);
...
// Add the n.Grad(u) boundary integral on the Dirichlet portion of the
// boundary marked in the dbc_marker array.
a.AddBdrFaceIntegrator(new DGDiffusionIntegrator(matCoef, sigma, kappa),
dbc_marker);
Where sigma
and kappa
are parameters controlling the symmetry and
interior penalty used in the DG diffusion formulation. These two
integrators work together to balance the natural boundary condition
associated with the DiffusionIntegrator
and to penalize solutions
which differ from the desired Dirichlet value near the boundary.
Similar pairs of integrators can be implemented to accommodate
other PDEs.