Nitsche method

this documentation is in english and needs to be translated.

1. Weak treatment of Dirichlet boundary conditions

In order to treat the boundary conditions uniformly (i.e. the same way as Neumann and Robin Conditions), we wish to treat the Dirichlet BC (e.g. $u=g$ ) weakly.

Initial Idea

Add the penalisation term

\int_{\partial Ω} μ (u - g)

$\int_{\partial \Omega} \mu( u - g)$

where $\mu$ is a constant. But this is not enough: this is not consistent with the initial formulation.

One can use the Nitsche method to implement weak Dirichlet conditions and follows the next steps:

write the equations in conservative form (i.e. identify the flux);
add the terms to ensure consistency (i.e the flux on the boundary);
symmetrize to ensure adjoint consistency;
add a penalisation term with factor $\gamma (u-g)/h$ that ensures that the solution will be set to the proper value at the boundary;

2. Penalisation parameter

2.1. Choosing $\gamma$

$\gamma$ must be chosen such that the coercivity(or inf-sup) property is satisfied. It is difficult to do in general. We increase $\gamma$ until the weak Dirichlet boundary conditions are properly satisfied, e.g. start with $\gamma=1$ , typical values are between 1 and 10.

The choice of $\gamma$ is a problem specially when $h$ is small.

2.2. Advantages, disadvantages

We compare the advantages and disadvantages of strong and weak Dirichlet boundary conditions treatment.

We start with the weak conditions

advantage uniform(weak) treatment of all boundary conditions type
advantage if boundary condition is independant of time, the terms are assembled once for all
disadvantage Introduction of the penalisation parameter $\gamma$ that needs to be tweaked

Strong treatment: Advantages

advantage Enforce exactely the boundary conditions
disadvantage Need to modify the matrix once assembled to reflect that the Dirichlet degree of freedom are actually known. Then even if the boundary condition is independant of time, at every time step if there are terms depending on time that need reassembly (e.g. convection) the strong treatment needs to be reapplied.
disadvantage it can be expensive to apply depending on the type of sparse matrix used, for example using CSR format setting rows to 0 except on the diagonal to 1 is not expensive but one must do that also for the columns associated with each Dirichlet degree of freedom and that is expensive.

3. Example: Laplacian

We look for $u$ such that

- \nabla \cdot (k \nabla u) = f in Ω, u = g |_{\partial Ω}

$-\nabla\cdot( k \nabla u )= f\ \mbox{in} \Omega,\quad u=g|_{\partial \Omega}$

$\begin{gathered} \label{eq:51} \int_\Omega k \nabla u \cdot \nabla v + \int_{\partial \Omega} \underbrace{-k \frac{\partial u}{\partial n}v}_{\text{integration by part}} \underbrace{-k \frac{\partial v}{\partial n} u}_{\text{adjoint consistency: symetrisation}} \\ + \underbrace{\frac{\gamma}{h} u v}_{\text{penalisation: enforce Dirichlet condition}} =\\ \int_\Omega f v + \int_{\partial \Omega} (\underbrace{-k \frac{\partial v}{\partial n} }_{\text{adjoint consistency}} + \underbrace{\frac{\gamma}{h} v) g}_{\text{penalisation: enforce Dirichlet condition}} \end{gathered}$

3.1. Implementation

// contribution to bilinear form (left hand side)
form2( _test=Xh, _trial=Xh ) +=
integrate( boundaryfaces(mesh),
           // integration by part
           -(gradt(u)*N())*id(v)
           // adjoint consistency
           -(grad(v)*N())*idt(u)
           // penalisation
           +gamma*id(v)*idt(u)/hFace());
// contribution to linear form (right hand side)
form1( _test=Xh ) +=
integrate( boundaryfaces(mesh),
           // adjoint consistency
           -(grad(v)*N())*g
           // penalisation
           +gamma*id(v)*g/hFace());

cpp

4. Example: Convection-Diffusion

Convection Diffusion Consider now the following problem, find $u$ such that

$\nabla \cdot ( -\epsilon\nabla u + \mathbf{c} u ) = f,\quad u = g|_{\partial \Omega},\quad \nabla \cdot \mathbf{c} = 0$

the flux vector field is $\mathbf{F}=-\nabla u$ .

Note that here the condition, $\nabla \cdot \mathbf{c} = 0$ was crucial to expand $\nabla \cdot (\mathbf{c} u )$ into $\mathbf{c} \cdot \nabla u$ since

$\nabla \cdot (\mathbf{c} u ) = \mathbf{c} \cdot \nabla u + \underbrace{u \nabla \cdot \mathbf{c}}_{=0}$

Weak formulation for convection diffusion Multiplying by any test function $v$ and integration by part of ([eq:2]) gives

$\int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} (\mathbf{F}\cdot \mathbf{n}) v = \int_\Omega f v$

where $\mathbf{n}$ is the outward unit normal to $\partial \Omega$ .

We now introduce the penalisation term that will ensure that $u \rightarrow g$ as $h \rightarrow 0$ on $\partial \Omega$ . ([eq:4]) reads now

$\int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} (\mathbf{F}\cdot \mathbf{n}) v + {\frac{\gamma}{h} u v} = \int_\Omega f v + {\int_{\partial \Omega} \frac{\gamma}{h} g v}$

Finally we can incorporate the symetrisation

$\begin{gathered} \int_\Omega \epsilon \nabla u \cdot \nabla v + (\mathbf{c} \cdot \nabla u)v + \int_{\partial \Omega} ((-\epsilon \nabla u)\cdot \mathbf{n}) v+ {(-\epsilon\nabla v\cdot \mathbf{n} + {\mathbf{c}}\cdot \mathbf{n} v ) u} + \frac{\gamma}{h} u v = \\ \int_\Omega f v + \int_{\partial \Omega} {(-\epsilon\nabla v\cdot \mathbf{n} + {\mathbf{c}}\cdot \mathbf{n} v) g}+ \frac{\gamma}{h} g v \end{gathered}$

4.1. Implementation

// bilinear form (left hand side)
form2( _trial=Xh, _test=Xh ) +=
integrate( boundaryfaces(mesh),
  // integration by part
  -($\epsilon$ gradt(u)*N())*id(v) + (idt(u)*trans(idv(c))*N())*id(v)
  // adjoint consistency
  -($\epsilon$ grad(v)*N())*idt(u) + (idt(u)*trans(idv(c))*N())*id(v)
  // penalisation
  +gamma*id(v)*idt(u)/hFace());
// linear form (right hand side)
form1( _test=Xh ) +=
integrate( boundaryfaces(mesh),
  // adjoint consistency
  -($\epsilon$ grad(v)*N())*g
  + g*trans(idv(c))*N())*id(v)
  // penalisation
  +gamma*id(v)*g/hFace());

cpp