Weak imposition of Dirichlet conditions for the Poisson problem#
Author: Jørgen S. Dokken
In this section, we will go through how to solve the Poisson problem from the Fundamentals tutorial using Nitsche’s method [Nit71]. The idea of weak imposition is that we add additional terms to the variational formulation to impose the boundary condition, instead of modifying the matrix system using strong imposition (lifting).
We start by importing the required modules and creating the mesh and function space for our solution
from dolfinx import fem, mesh, plot, default_scalar_type
from dolfinx.fem.petsc import LinearProblem
import numpy
from mpi4py import MPI
from ufl import (Circumradius, FacetNormal, SpatialCoordinate, TrialFunction, TestFunction,
div, dx, ds, grad, inner)
N = 8
domain = mesh.create_unit_square(MPI.COMM_WORLD, N, N)
V = fem.functionspace(domain, ("Lagrange", 1))
Next, we create a function containing the exact solution (which will also be used in the Dirichlet boundary condition) and the corresponding source function for the right hand side. Note that we use ufl.SpatialCoordinate
to define the exact solution, which in turn is interpolated into uD
and used to create the source function f
.
uD = fem.Function(V)
x = SpatialCoordinate(domain)
u_ex = 1 + x[0]**2 + 2 * x[1]**2
uD.interpolate(fem.Expression(u_ex, V.element.interpolation_points()))
f = -div(grad(u_ex))
As opposed to the first tutorial, we now have to have another look at the variational form. We start by integrating the problem by parts, to obtain
As we are not using strong enforcement, we do not set the trace of the test function to \(0\) on the outer boundary. Instead, we add the following two terms to the variational formulation
where the first term enforces symmetry to the bilinear form, while the latter term enforces coercivity. \(u_D\) is the known Dirichlet condition, and \(h\) is the diameter of the circumscribed sphere of the mesh element. We create bilinear and linear form, \(a\) and \(L\)
u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(domain)
h = 2 * Circumradius(domain)
alpha = fem.Constant(domain, default_scalar_type(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
As we now have the variational form, we can solve the linear problem
problem = LinearProblem(a, L)
uh = problem.solve()
We compute the error of the computation by comparing it to the analytical solution
error_form = fem.form(inner(uh-uD, uh-uD) * dx)
error_local = fem.assemble_scalar(error_form)
errorL2 = numpy.sqrt(domain.comm.allreduce(error_local, op=MPI.SUM))
if domain.comm.rank == 0:
print(fr