对于 H(div) 有限元,物理单元上的向量场不能随意从参考元“搬过来”。
设参考元到物理元的映射为
F:K^→K,x=F(x^),
其 Jacobian 为
J=DF(x^).
在 H(div) 元中,正确的变换不是普通的函数复合,而是 Piola 变换:
u(x)=detJ1Ju^(x^),x=F(x^).
其中 u^ 是参考元上的基函数。
问题是:为什么一定要这样定义?为什么不能直接把参考元基函数“拉回来”定义为
u(x)=u^(F−1(x))?
答案在于:H(div) 元需要保持的是法向通量和散度结构,而 Piola 变换正是保持这两个量的自然变换。
与 H1 元不同,H(div) 元最重要的不是点值,而是:
边上的法向通量自由度
∫eu⋅nqds.
对 RT0 元来说,q 只是常数,所以本质上就是每条边上的总法向通量。
散度的正确变换
∇⋅u.
因此,物理元上的基函数必须满足:
- 参考元上的边通量自由度映到物理元后仍然保持一致;
- 参考元上的常散度映到物理元后仍然是正确的常散度;
- 参考元上的 RT 空间映到物理元后仍然落在物理元上的 RT 空间中。
Piola 变换正好满足这些要求。
对仿射映射 x=F(x^)=Jx^+b,Piola 变换有两个核心性质。
若 e=F(e^),则有
∫eu⋅nds=∫e^u^⋅n^ds^.
这说明边通量自由度在映射后严格保持。
∇x⋅u=detJ1∇x^⋅u^.
因此参考元上常散度的函数,映到物理元后仍然是常散度函数。
取参考三角形
K^=conv{(0,0),(1,0),(0,1)}.
它的三条边分别是:
- e^1:η=0,底边;
- e^2:ξ+η=1,斜边;
- e^3:ξ=0,左边。
取 RT0 的三个基函数为
ϕ^1=(ξη−1),ϕ^2=(ξη),ϕ^3=(ξ−1η).
它们分别对应三条边的法向通量自由度。
在底边 e^1:η=0 上,外法向可取
n^1=(0,−1).
此时
ϕ^1⋅n^1=(ξ,−1)⋅(0,−1)=1.
所以它在底边上的法向通量为 1。
在左边 e^3:ξ=0 上,外法向可取
n^3=(−1,0),
则
ϕ^1⋅n^3=(0,η−1)⋅(−1,0)=0.
在斜边 e^2:ξ+η=1 上,外法向可取
n^2=(1,1),
于是
ϕ^1⋅n^2=ξ+(η−1)=ξ+η−1=0.
所以 ϕ^1 只在第一条边上有单位通量。
类似地:
- ϕ^2 对应斜边;
- ϕ^3 对应左边。
取仿射映射
F(ξ,η)=(xy)=(2011)(ξη).
也就是
x=2ξ+η,y=η.
于是物理三角形 K 的顶点为
(0,0),(2,0),(1,1).
Jacobian 为
J=(2011),detJ=2.
逆映射为
η=y,ξ=2x−y.
定义
ϕi(x,y)=21Jϕ^i(ξ,η),(ξ,η)=F−1(x,y).
ϕ^1=(ξη−1).
先乘以 J:
Jϕ^1=(2011)(ξη−1)=(2ξ+η−1η−1).
再除以 detJ=2:
ϕ1=21(2ξ+η−1η−1).
由于 x=2ξ+η, y=η,可写成
ϕ1(x,y)=21(x−1y−1)
ϕ^2=(ξη)⟹ϕ2=21(2ξ+ηη)=21(xy)
ϕ^3=(ξ−1η)⟹ϕ3=21(2ξ+η−2η)=21(x−2y)
物理三角形上的 RT0 空间函数具有形式
v(x,y)=(ab)+c(xy).
而我们得到的三个函数为
\boldsymbol\phi_1= rac12\binom{x-1}{y-1}, \qquad \boldsymbol\phi_2= rac12\binom{x}{y}, \qquad \boldsymbol\phi_3= rac12\binom{x-2}{y}.
它们都恰好是
(a,b)+c(x,y)
这种形式,因此确实仍然属于物理元上的 RT0 空间。
这条边满足 y=0,外法向取
n1=(0,−1).
对 ϕ1,有
ϕ1(x,0)=21(x−1−1).
因此
ϕ1⋅n1=21.
由于底边长度为 2,所以总通量为
∫e1ϕ1⋅n1ds=21×2=1.
恰好保持参考元上的自由度。
同时
ϕ2⋅n1=0,ϕ3⋅n1=0.
所以第一条边的通量自由度结构完全正确。
这条边方程是 x+y=2,可取外法向
n2=(1,1).
对 ϕ2,有
ϕ2⋅n2=21(x+y).
在边 e2 上,x+y=2,故
ϕ2⋅n2=1.
因此它在第二条边上有单位通量,另外两个基函数在此边上法向通量为 0。
这条边方程是 x=y,可取外法向
n3=(−1,1).
对 ϕ3,有
ϕ3⋅n3=21(−(x−2)+y)=21(2−x+y).
在 x=y 上有
ϕ3⋅n3=1.
所以第三个基函数在第三条边上具有单位通量。
综上,Piola 变换后的三个基函数仍然满足 RT0 元“每条边一个基函数对应单位通量,其余边通量为 0”的定义。
现在不用 Piola,而直接定义
ϕ~i(x,y)=ϕ^i(F−1(x,y)).
以第一基函数为例:
ϕ~1(x,y)=(ξη−1)=(2x−yy−1).
RT0 函数应具有形式
(a,b)+c(x,y).
但这里
ϕ~1(x,y)=(21x−21yy−1).
第一分量里既有 x 又有 y,第二分量中又是另一个不同系数的 y,这不可能写成
(a,b)+c(x,y)
的形式,所以它已经不在物理单元的 RT0 空间里。
在底边 e1:y=0 上,
ϕ~1(x,0)=(2x−1).
用外法向 n1=(0,−1),可得
ϕ~1⋅n1=1.
注意:底边长度现在是 2,因此总通量为
∫e1ϕ~1⋅n1ds=1×2=2.
而参考元对应的自由度是 1。也就是说,直接拉回以后,自由度被放大了一倍。
所以它根本不能作为正确的 H(div) 基函数。
参考元上
∇x^⋅ϕ^1=∂ξ∂ξ+∂η∂(η−1)=2.
Piola 变换后,理论上应满足
∇x⋅ϕ1=detJ1⋅2=1.
而我们确实有
ϕ1=21(y−1x−1),
故
∇x⋅ϕ1=21+21=1.
完全正确。
但如果直接拉回,
ϕ~1=(y−121x−21y),
则
∇x⋅ϕ~1=21+1=23,
这已经不是正确的变换结果了。
因此,H(div) 有限元必须使用 Piola 变换,而不能直接把参考元基函数做普通拉回。原因是:
- Piola 变换保持边法向通量自由度不变;
- Piola 变换保持散度的正确变换规律;
- Piola 变换把参考元上的 RT 空间映成物理元上的 RT 空间;
- 只有这样,单元间法向连续性才有意义,H(div) 元的定义才成立。
而普通拉回会导致:
- 法向通量自由度错误;
- 散度错误;
- 映射后的函数不再属于物理元上的 RT0 空间。
所以对于 H(div) 元,Piola 变换不是一种“方便的选法”,而是由有限元结构本身逼出来的唯一正确选择。
H(div) 元要保持的是“法向通量”,不是“向量值本身”;Piola 变换正是保持法向通量和散度结构的变换,因此 H(div) 基函数必须通过 Piola 变换定义。