Weak imposition of Dirichlet conditions for the Poisson problem

Weak imposition of Dirichlet conditions for the Poisson problem#

Author: Jørgen S. Dokken 翻译：

章节使用Nitsche方法 [Nit71]求解章节Fundamentals中的Poisson问题。 (Weak imposition是什么意思？ strong imposition：修改矩阵元素，也称为lifting。 weak imposition：通过给变分公式添加额外的项来施加边界条件。）

首先，我们定义网格和函数空间。

from dolfinx import fem, mesh, plot
import numpy
from mpi4py import MPI
from petsc4py.PETSc import ScalarType 
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, ("CG", 1))

接下来，我们定义一个真实解和对应的源项。用 ufl.SpatialCoordinate 定义真实解，对它插值得到边界条件 uD ，并推导出源项 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 与第一章的教程不同，我们变分形式。得到

(8)#\[\begin{align} \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d}x - \int_{\partial\Omega}\nabla u \cdot n v~\mathrm{d}s = \int_{\Omega} f v~\mathrm{d}x. \end{align}\]

因为我们没有施加强制边界条件，所以我们没有在外部边界上将测试函数的迹设为 \(0\) 。取而代之，我们在变分公式中加入下面两项。

(9)#\[\begin{align} -\int_{\partial\Omega} \nabla v \cdot n (u-u_D)~\mathrm{d}s + \frac{\alpha}{h} \int_{\partial\Omega} (u-u_D)v~\mathrm{d}s. \end{align}\]

其中，第一项强制双线性型的对称性，第二项为矫顽力（ coercivity ?）。\(u_D\) 为Dirichlet边界条件， \(h\) 为网格单元外接球的直径。双线性型 \(a\) 和线性型 \(L\)分别为

(10)#\[\begin{align} a(u, v) &= \int_{\Omega} \nabla u \cdot \nabla v~\mathrm{d}x + \int_{\partial\Omega}-(n \cdot\nabla u) v - (n \cdot \nabla v) u + \frac{\alpha}{h} uv~\mathrm{d}s,\\ L(v) &= \int_{\Omega} fv~\mathrm{d}x + \int_{\partial\Omega} -(n \cdot \nabla v) u_D + \frac{\alpha}{h} u_Dv~\mathrm{d}s \end{align}\]

u = TrialFunction(V)
v = TestFunction(V)
n = FacetNormal(domain)
h = 2 * Circumradius(domain)
alpha = fem.Constant(domain, ScalarType(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 = fem.petsc.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"$L^2$-error: {errorL2:.2e}")

$L^2$-error: 1.59e-03

We observe that the \(L^2\)-error is of the same magnitude as in the first tutorial. As in the previous tutorial, we also compute the maximal error for all the degrees of freedom.

error_max = domain.comm.allreduce(numpy.max(numpy.abs(uD.x.array-uh.x.array)), op=MPI.MAX)
if domain.comm.rank == 0:
    print(f"Error_max : {error_max:.2e}")

Error_max : 5.41e-03

We observe that as we weakly impose the boundary condition, we no longer fullfill the equation to machine precision at the mesh vertices. We also plot the solution using pyvista

import pyvista
pyvista.start_xvfb()

grid = pyvista.UnstructuredGrid(*plot.create_vtk_mesh(V))
grid.point_data["u"] = uh.x.array.real
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_mesh(grid, show_edges=True, show_scalar_bar=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
    plotter.show()
else:
    figure = plotter.screenshot("nitsche.png")

Nit71: J. Nitsche. Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, 36(1):9–15, Jul 1971. doi:10.1007/BF02995904.