Nitsche's method for imposition of Dirichlet boundary conditions: implementation standpoint

I'm trying to understand how Nitsche's method works in practice. I understood the theoretical principle behind it, but what I can't understand is its implementation. More precisely, I'd like to solve the classical Poisson equation on the square with classical conforming degree 1 finite elements using Nitsche's method for the imposition of boundary conditions.

The variational formulation is to find $$u_h \in V_h$$ s.t.

$$(\nabla u_h, \nabla v ) - \langle \nabla u_h \cdot n, v \rangle-\langle u,\nabla v \cdot n \rangle + \sum_{E \in \text{ boundary faces}}\beta h_E^{-1} $$ $$=$$ $$(f,v)- \langle u_0, \nabla v\cdot n\rangle + \sum_{E \in \text{ boundary faces}} \beta h_E^{-1} \langle u_0,v \rangle$$

for every $$v \in V_h$$. The formulation is taken from here, pag.2.

As you can see, all the $$L^2$$ inner products are taken on the boundary faces, no interior ones.

My guess for the implementation standpoint is: one doesn't have to post-process the matrix of the linear system, i.e. one has to stop after all the local contributions have been distributed into the global matrix. Indeed, looking at the formulation, what is different from the usual $$(\nabla u, \nabla v)=(f,v)$$ weak form are all terms on the boundary. Is that correct?

You are correct, all degrees-of-freedom are constrained weakly so there is no need to post process the matrix. Here is an example with $$f=10$$ and $$u_0(x) = \sin(2 \pi x)$$: The example source code runs after pip install scikit-fem==4.0.1:

import numpy as np
from skfem import *
from skfem.visuals.matplotlib import plot, show

m = MeshTri.init_sqsymmetric().refined(4)
e = ElementTriP1()
alpha = 1e-3

ib = Basis(m, e)
bb = FacetBasis(m, e)

def u0(x):
return np.sin(2. * np.pi * x)

@BilinearForm
def nitsche_bilinf(u, v, p):
h = p.h
n = p.n
return u * v / (alpha * h) - dot(grad(u), n) * v - dot(grad(v), n) * u

@LinearForm
h = p.h
n = p.n
return u0(p.x) * v / (alpha * h) - u0(p.x) * dot(grad(v), n)

A = asm(laplace, ib)
B = asm(nitsche_bilinf, bb)
f = 10 * asm(unit_load, ib)

x = solve(A + B, f + g)

plot(ib, x, colorbar=True)
show()

You can experiment with the source code also in Google Colab.

• Thanks for your answer, much appreciated. May I ask if the solution you obtain is the same as the one you have by post-processing the matrix? As a matter of fact, if I test my method with $\beta = 10.0$ and homogeneous Dirichlet and $f=1$, I have as boundary values something like $10^{-5}$, which seems not satisfactory at all Nov 8 '21 at 9:46
• The solution is not the same. However, the convergence rate towards the exact solution should be the same. In theory, there exists a specific optimal value for $\beta$ which can be defined as a solution to a specific eigenvalue problem. However, this is a topic for another time.
– knl
Nov 8 '21 at 14:49
• Let me add something that you probably already noticed: if you increase $\beta$ (in my formulation $\beta = \alpha^{-1}$ so you need to decrease $\alpha$) then your solution on the boundary is closer to zero in the homogeneous case.
– knl
Nov 8 '21 at 14:55
• That's right. Indeed, I'm a bit puzzled by the fact that with that test case, using $66049$ DoFs my solution has boundary values of the order $10^{-7}$ and and the maximum value is ca. $7.4 \cdot 10^{-2}$, instead of $0.29$, as it should be. Do you think there are some problems in the implementation, or is this reasonable? (I'm using deal.ii library). @knl Nov 8 '21 at 15:45
• Okay, I've just verified that indeed the values are really pretty much equal in the case $f=1$ with homogeneous Dirichlet by using your colab page (btw, thanks!) Just one thing: how is it possible that your minimum is about $-2e-6$ and mine $+2e-6$? That's strange, may it depend on the solver your library is using? Nov 8 '21 at 16:32

If I understand your question right then yes, you're correct. The most common approach to enforcing Dirichlet boundary conditions with the finite element method is to modify the linear system of equations, which could be called a post-processing step after matrix assembly to paraphrase you. Nitsche's method circumvents the need for this post-processing step by instead modifying the variational principle. The key difference with Nitsche's method is that now there are some integrals over the boundary in the weak form that must also be assembled into the matrix. In terms of complexity of implementation, adding these extra boundary integrals is about the same difficulty as adding Robin boundary conditions.

One of the advantages of Nitsche is that the way that you modify the problem to be solved is pretty much universal across implementations, whereas the necessary post-processing of the matrix is going to depend on the individual software package you use for sparse linear algebra. If you want to read more, I wrote a bit about the motivation for Nitsche's method and how you can derive a lower bound for the penalty parameter $$\beta$$ here.

Another advantage is that, while some boundary conditions can be implemented by tweaking a few matrix and RHS vector entries, other types of BCs cannot. For example, Stokes flow with frictional slip (rather than no slip) is part Dirichlet, part Robin boundary conditions. If the boundaries aren't a straight line, the usual approach won't work at all, but Nitsche's method is quite easy.

• I think your first blog post is an interesting overview but may not be exactly accurate on the Lagrange multiplier approach and the contribution of Pitkäranta, e.g., based on my understanding of the papers and a quick test P1-P1 may work just fine in the usual cases. Nitsche's method is exactly the residual stabilization applied to the Lagrange multiplier approach (see Stenberg - On some techniques for approximating boundary conditions in the finite element method) which is also a useful classification leading to rigorous numerical analysis.
– knl
Nov 7 '21 at 19:16
• @knl I always thought that using Lagrange multipliers for essential BCs just wasn't done very much because Pitkäranta showed that satisfying LBB for the multiplier space is so difficult as to be infeasible. Did I misunderstand Pitkäranta's results? Or did I read what he said correctly but you think some cases can actually work out fine? Nov 8 '21 at 18:07
• Here is an implementation of P1-P1 method using the trace mesh for the Lagrange multiplier: colab.research.google.com/drive/… You can try it out for different meshes and boundary data. It doesn't tell anything about uniform stability but I wasn't able to find any obvious cases where the matrix would be singular. Uniform stability is another story.
– knl
Nov 10 '21 at 11:10
• Now reading about description of the issue in Barbosa-Hughes (doi:10.1016/0045-7825(91)90125-p) confuse me because they say equal order interpolation does not work. Maybe I still didn't do a complex enough example?
– knl
Nov 10 '21 at 11:13
• I note that another application of Nitsche's method (or in general the weak imposition of BCs) is non-conforming meshes. In that case, you cannot satisfy the Dirichlet BC by imposing the values at the nodes. Nov 13 '21 at 21:55