Nitsche method for dirichlet boundary conditions

variational formulation
1

Hello. I am attempting to implement Nitsche method for general boundary conditions from a paper by Juntunen and Stenberg, and my first step is to reproduce Nitsche method for Dirichlet boundary conditions. I found an old tutorial by Dokkens and have successfully translated the code to latest version of dolfinx. However, I’m struggling to understand the origin of the term u_D and the linked “fundamentals” markdown is broken so I cannot see the definition of the problem that is being solved in the nitsche notebook. Could someone guide me to an implementation of Nitsche’s method for simple Dirichlet boundary conditions (e.g. u(x=0, y) = 1, u(x=Lx, y) = 0 and u(x, y=0) = u(x, y=Ly) = 0)?

Here’s my translation of the Nitsche’s demo in the linked notebook but without the pyvista visualization.

import dolfinx
import numpy as np
import pyvista
import ufl

from dolfinx import cpp, fem, io, mesh, plot
from mpi4py import MPI
from petsc4py import PETSc
from ufl import (div, dx, ds, grad, inner, grad)

comm = MPI.COMM_WORLD
N = 8
domain = mesh.create_unit_square(comm, N, N)
V = fem.FunctionSpace(domain, ("CG", 1))

uD = fem.Function(V)
uD.interpolate(lambda x: 1 + x[0] ** 2 + 2 * x[1] ** 2)
uD.x.scatter_forward()
x = ufl.SpatialCoordinate(domain)
n = ufl.FacetNormal(domain)

f = -div(grad(1 + x[0] ** 2 + 2 * x[1] ** 2))

u = ufl.TrialFunction(V)
v = ufl.TestFunction(V)
h = 2 * ufl.Circumradius(domain)
alpha = 10
a = inner(grad(u), grad(v)) * dx - inner(n, grad(u)) * v * ds
a += - inner(n, grad(v)) * u * ds + alpha / h * inner(u, v) * ds
L = inner(f, v) * dx
L += - inner(n, grad(v)) * uD * ds + alpha / h * inner(uD, v) * ds

problem = fem.petsc.LinearProblem(a, L)
uh = problem.solve()

with io.XDMFFile(comm, "nitsche.xdmf", "w") as outfile:
    outfile.write_mesh(domain)
    outfile.write_function(uh)

error_form = fem.form(inner(uh - uD, uh - uD) * dx)
errorL2 = np.sqrt(fem.assemble_scalar(error_form))
print(fr"$L^2$-error: {errorL2:.2e}")

u_vertex_values = uh.x.array.real
u_ex_vertex_values = uD.x.array.real
error_max = np.max(np.abs(u_vertex_values - u_ex_vertex_values))
print(f"Error_max : {error_max:.2e}")

Thank you for any help.

2

The fundamentals demon in question is simply:
https://jsdokken.com/dolfinx-tutorial/chapter1/fundamentals.html

As it is solving a problem with u=u_D on the boundary, you need to add a function with these values (or use ufl to define the bc, either with constants or spatial-coordinates).

You can either spilt the boundary conditions into four similar integrals over each sub component of the boundary. See for instance Setting multiple Dirichlet, Neumann, and Robin conditions — FEniCSx tutorial

3

Thank you! That worked.