Nitsche 施加狄利克雷边界条件的方法:实施观点

计算科学 有限元 交易.ii 弱解 尼采法
2021-11-28 18:53:43

我试图了解 Nitsche 的方法在实践中是如何工作的。我理解它背后的理论原理,但我无法理解的是它的实现。更准确地说,我想使用 Nitsche 施加边界条件的方法来求解具有经典符合度 1 有限元的正方形上的经典 Poisson 方程。

变分公式是找到uhVh英石

(uh,v)uhn,vu,vn+E boundary facesβhE1<uh,v>
=
(f,v)u0,vn+E boundary facesβhE1u0,v

对于每个vVh. 该公式取自此处,第 2 页

如您所见,所有L2内积取在边界面上,没有内积。

我对实现观点的猜测是:不必对线性系统的矩阵进行后处理,即必须在所有局部贡献都分配到全局矩阵后停止。确实,看配方,与通常的不同(u,v)=(f,v)弱形式都是边界上的项。那是对的吗?

2个回答

如果我理解你的问题,那么是的,你是对的。使用有限元方法强制执行狄利克雷边界条件的最常见方法是修改线性方程组,这可以称为矩阵组装后的后处理步骤来解释您的意思。 Nitsche 的方法通过修改变分原理来规避此后处理步骤的需要。 与 Nitsche 方法的主要区别在于,现在边界上有一些弱形式的积分也必须组装到矩阵中。就实现的复杂性而言,添加这些额外的边界积分与添加 Robin 边界条件的难度大致相同。

Nitsche 的优点之一是修改要解决的问题的方式在实现中几乎是通用的,而矩阵的必要后处理将取决于您用于稀疏线性代数的单个软件包。如果您想阅读更多内容,我写了一些关于 Nitsche 方法的动机以及如何得出惩罚参数的下限β 在这里

另一个优点是,虽然可以通过调整一些矩阵和 RHS 向量条目来实现一些边界条件,但其他类型的 BC 不能。例如,具有摩擦滑移(而不是无滑移)的斯托克斯流是部分狄利克雷,部分罗宾边界条件。如果边界不是一条直线,通常的方法根本行不通,但 Nitsche 的方法很简单

你是对的,所有自由度都受到弱约束,因此无需对矩阵进行后处理。这是一个例子f=10u0(x)=sin(2πx)

数值解

示例源代码在以下之后运行pip install scikit-fem==4.0.1

import numpy as np
from skfem import *
from skfem.helpers import grad, dot
from skfem.models import laplace, unit_load
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[0])

@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
def nitsche_load(v, p):
    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)
g = asm(nitsche_load, bb)

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

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

您也可以在Google Colab中试验源代码