我的博客

外观

hdiv-piola

约 2283 字大约 8 分钟

2026-03-23

1. 问题背景

对于 H(div)H(\mathrm{div}) 有限元，物理单元上的向量场不能随意从参考元“搬过来”。

设参考元到物理元的映射为

F:K^→K,x=F(x^),F:\hat K\to K,\qquad x = F(\hat x),

其 Jacobian 为

J=DF(x^).J = DF(\hat x).

在 H(div)H(\mathrm{div}) 元中，正确的变换不是普通的函数复合，而是 Piola 变换：

u(x)=1det⁡JJu^(x^),x=F(x^).\mathbf u(x)=\frac{1}{\det J}J\hat{\mathbf u}(\hat x), \qquad x=F(\hat x).

其中 u^\hat{\mathbf u} 是参考元上的基函数。

问题是：为什么一定要这样定义？为什么不能直接把参考元基函数“拉回来”定义为

u(x)=u^(F−1(x))?\mathbf u(x)=\hat{\mathbf u}(F^{-1}(x))?

答案在于：H(div)H(\mathrm{div}) 元需要保持的是法向通量和散度结构，而 Piola 变换正是保持这两个量的自然变换。

2. H(div)H(\mathrm{div}) 元真正要保持的是什么

与 H1H^1 元不同，H(div)H(\mathrm{div}) 元最重要的不是点值，而是：

  1. 边上的法向通量自由度

    ∫eu⋅n q ds.\int_e \mathbf u\cdot \mathbf n\, q\, ds.

    对 RT0 元来说，qq 只是常数，所以本质上就是每条边上的总法向通量。

  2. 散度的正确变换

    ∇⋅u.\nabla\cdot \mathbf u.

因此，物理元上的基函数必须满足：

  • 参考元上的边通量自由度映到物理元后仍然保持一致；
  • 参考元上的常散度映到物理元后仍然是正确的常散度；
  • 参考元上的 RT 空间映到物理元后仍然落在物理元上的 RT 空间中。

Piola 变换正好满足这些要求。

3. Piola 变换的两个关键性质

对仿射映射 x=F(x^)=Jx^+bx=F(\hat x)=J\hat x+b，Piola 变换有两个核心性质。

3.1 保持边法向通量

若 e=F(e^)e=F(\hat e)，则有

∫eu⋅n ds=∫e^u^⋅n^ ds^.\int_e \mathbf u\cdot \mathbf n\, ds = \int_{\hat e}\hat{\mathbf u}\cdot \hat{\mathbf n}\, d\hat s.

这说明边通量自由度在映射后严格保持。

3.2 散度按正确公式变换

∇x⋅u=1det⁡J∇x^⋅u^.\nabla_x\cdot \mathbf u = \frac{1}{\det J}\nabla_{\hat x}\cdot \hat{\mathbf u}.

因此参考元上常散度的函数，映到物理元后仍然是常散度函数。

4. 参考三角形上的 RT0 元

取参考三角形

K^=conv{(0,0),(1,0),(0,1)}.\hat K = \mathrm{conv}\{(0,0),(1,0),(0,1)\}.

它的三条边分别是：

  • e^1:η=0\hat e_1: \eta = 0，底边；
  • e^2:ξ+η=1\hat e_2: \xi+\eta=1，斜边；
  • e^3:ξ=0\hat e_3: \xi = 0，左边。

取 RT0 的三个基函数为

ϕ^1=(ξη−1),ϕ^2=(ξη),ϕ^3=(ξ−1η).\hat{\boldsymbol\phi}_1= \begin{pmatrix} \xi\\ \eta-1 \end{pmatrix}, \qquad \hat{\boldsymbol\phi}_2= \begin{pmatrix} \xi\\ \eta \end{pmatrix}, \qquad \hat{\boldsymbol\phi}_3= \begin{pmatrix} \xi-1\\ \eta \end{pmatrix}.

它们分别对应三条边的法向通量自由度。

5. 参考元上法向通量的验证

5.1 第一基函数 ϕ^1\hat{\boldsymbol\phi}_1

在底边 e^1:η=0\hat e_1: \eta=0 上，外法向可取

n^1=(0,−1).\hat n_1=(0,-1).

此时

ϕ^1⋅n^1=(ξ,−1)⋅(0,−1)=1.\hat{\boldsymbol\phi}_1\cdot \hat n_1 =(\xi,-1)\cdot(0,-1)=1.

所以它在底边上的法向通量为 1。

在左边 e^3:ξ=0\hat e_3: \xi=0 上，外法向可取

n^3=(−1,0),\hat n_3=(-1,0),

则

ϕ^1⋅n^3=(0,η−1)⋅(−1,0)=0.\hat{\boldsymbol\phi}_1\cdot \hat n_3 =(0,\eta-1)\cdot(-1,0)=0.

在斜边 e^2:ξ+η=1\hat e_2: \xi+\eta=1 上，外法向可取

n^2=(1,1),\hat n_2=(1,1),

于是

ϕ^1⋅n^2=ξ+(η−1)=ξ+η−1=0.\hat{\boldsymbol\phi}_1\cdot \hat n_2 = \xi+(\eta-1)=\xi+\eta-1=0.

所以 ϕ^1\hat{\boldsymbol\phi}_1 只在第一条边上有单位通量。

类似地：

  • ϕ^2\hat{\boldsymbol\phi}_2 对应斜边；
  • ϕ^3\hat{\boldsymbol\phi}_3 对应左边。

6. 构造一个物理单元

取仿射映射

F(ξ,η)=(xy)=(2101)(ξη).F(\xi,\eta)= \begin{pmatrix} x\\y \end{pmatrix} = \begin{pmatrix} 2 & 1\\ 0 & 1 \end{pmatrix} \begin{pmatrix} \xi\\\eta \end{pmatrix}.

也就是

x=2ξ+η,y=η.x=2\xi+\eta, \qquad y=\eta.

于是物理三角形 KK 的顶点为

(0,0),(2,0),(1,1).(0,0),\qquad (2,0),\qquad (1,1).

Jacobian 为

J=(2101),det⁡J=2.J= \begin{pmatrix} 2 & 1\\ 0 & 1 \end{pmatrix}, \qquad \det J=2.

逆映射为

η=y,ξ=x−y2.\eta=y, \qquad \xi=\frac{x-y}{2}.

7. 用 Piola 变换构造物理元基函数

定义

ϕi(x,y)=12Jϕ^i(ξ,η),(ξ,η)=F−1(x,y).\boldsymbol\phi_i(x,y)=\frac{1}{2}J\hat{\boldsymbol\phi}_i(\xi,\eta), \qquad (\xi,\eta)=F^{-1}(x,y).

7.1 第一基函数

ϕ^1=(ξη−1).\hat{\boldsymbol\phi}_1= \begin{pmatrix} \xi\\ \eta-1 \end{pmatrix}.

先乘以 JJ：

Jϕ^1=(2101)(ξη−1)=(2ξ+η−1η−1).J\hat{\boldsymbol\phi}_1= \begin{pmatrix} 2 & 1\\ 0 & 1 \end{pmatrix} \begin{pmatrix} \xi\\ \eta-1 \end{pmatrix} = \begin{pmatrix} 2\xi+\eta-1\\ \eta-1 \end{pmatrix}.

再除以 det⁡J=2\det J=2：

ϕ1=12(2ξ+η−1η−1).\boldsymbol\phi_1 = \frac12 \begin{pmatrix} 2\xi+\eta-1\\ \eta-1 \end{pmatrix}.

由于 x=2ξ+ηx=2\xi+\eta, y=ηy=\eta，可写成

ϕ1(x,y)=12(x−1y−1)\boxed{ \boldsymbol\phi_1(x,y)= \frac12 \begin{pmatrix} x-1\\ y-1 \end{pmatrix}}

7.2 第二基函数

ϕ^2=(ξη)⟹ϕ2=12(2ξ+ηη)=12(xy)\hat{\boldsymbol\phi}_2= \begin{pmatrix} \xi\\ \eta \end{pmatrix} \quad\Longrightarrow\quad \boldsymbol\phi_2 = \frac12 \begin{pmatrix} 2\xi+\eta\\ \eta \end{pmatrix} = \boxed{ \frac12 \begin{pmatrix} x\\y \end{pmatrix}}

7.3 第三基函数

ϕ^3=(ξ−1η)⟹ϕ3=12(2ξ+η−2η)=12(x−2y)\hat{\boldsymbol\phi}_3= \begin{pmatrix} \xi-1\\ \eta \end{pmatrix} \quad\Longrightarrow\quad \boldsymbol\phi_3 = \frac12 \begin{pmatrix} 2\xi+\eta-2\\ \eta \end{pmatrix} = \boxed{ \frac12 \begin{pmatrix} x-2\\y \end{pmatrix}}

8. Piola 变换后的函数仍然属于 RT0 空间

物理三角形上的 RT0 空间函数具有形式

v(x,y)=(ab)+c(xy).\mathbf v(x,y)= \begin{pmatrix} a\\b \end{pmatrix} + c \begin{pmatrix} x\\y \end{pmatrix}.

而我们得到的三个函数为

\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)(a,b)+c(x,y)

这种形式，因此确实仍然属于物理元上的 RT0 空间。

9. 验证 Piola 基函数的法向通量

9.1 边 e1e_1：底边 (0,0)→(2,0)(0,0)\to(2,0)

这条边满足 y=0y=0，外法向取

n1=(0,−1).n_1=(0,-1).

对 ϕ1\boldsymbol\phi_1，有

ϕ1(x,0)=12(x−1−1).\boldsymbol\phi_1(x,0) = \frac12 \begin{pmatrix} x-1\\ -1 \end{pmatrix}.

因此

ϕ1⋅n1=12.\boldsymbol\phi_1\cdot n_1 = \frac12.

由于底边长度为 2，所以总通量为

∫e1ϕ1⋅n1 ds=12×2=1.\int_{e_1}\boldsymbol\phi_1\cdot n_1\,ds = \frac12\times 2=1.

恰好保持参考元上的自由度。

同时

ϕ2⋅n1=0,ϕ3⋅n1=0.\boldsymbol\phi_2\cdot n_1=0, \qquad \boldsymbol\phi_3\cdot n_1=0.

所以第一条边的通量自由度结构完全正确。

9.2 边 e2e_2：斜边 (2,0)→(1,1)(2,0)\to(1,1)

这条边方程是 x+y=2x+y=2，可取外法向

n2=(1,1).n_2=(1,1).

对 ϕ2\boldsymbol\phi_2，有

ϕ2⋅n2=12(x+y).\boldsymbol\phi_2\cdot n_2 = \frac12(x+y).

在边 e2e_2 上，x+y=2x+y=2，故

ϕ2⋅n2=1.\boldsymbol\phi_2\cdot n_2=1.

因此它在第二条边上有单位通量，另外两个基函数在此边上法向通量为 0。

9.3 边 e3e_3：左边 (0,0)→(1,1)(0,0)\to(1,1)

这条边方程是 x=yx=y，可取外法向

n3=(−1,1).n_3=(-1,1).

对 ϕ3\boldsymbol\phi_3，有

ϕ3⋅n3=12(−(x−2)+y)=12(2−x+y).\boldsymbol\phi_3\cdot n_3 = \frac12\bigl(-(x-2)+y\bigr) = \frac12(2-x+y).

在 x=yx=y 上有

ϕ3⋅n3=1.\boldsymbol\phi_3\cdot n_3=1.

所以第三个基函数在第三条边上具有单位通量。

综上，Piola 变换后的三个基函数仍然满足 RT0 元“每条边一个基函数对应单位通量，其余边通量为 0”的定义。

10. 如果直接拉回，会发生什么

现在不用 Piola，而直接定义

ϕ~i(x,y)=ϕ^i(F−1(x,y)).\tilde{\boldsymbol\phi}_i(x,y)=\hat{\boldsymbol\phi}_i(F^{-1}(x,y)).

以第一基函数为例：

ϕ~1(x,y)=(ξη−1)=(x−y2y−1).\tilde{\boldsymbol\phi}_1(x,y) = \begin{pmatrix} \xi\\ \eta-1 \end{pmatrix} = \begin{pmatrix} \dfrac{x-y}{2}\\ y-1 \end{pmatrix}.

10.1 它已经不属于物理元上的 RT0 空间

RT0 函数应具有形式

(a,b)+c(x,y).(a,b)+c(x,y).

但这里

ϕ~1(x,y)=(12x−12yy−1).\tilde{\boldsymbol\phi}_1(x,y) = \begin{pmatrix} \frac12 x-\frac12 y\\ y-1 \end{pmatrix}.

第一分量里既有 xx 又有 yy，第二分量中又是另一个不同系数的 yy，这不可能写成

(a,b)+c(x,y)(a,b)+c(x,y)

的形式，所以它已经不在物理单元的 RT0 空间里。

10.2 它的边通量也不对

在底边 e1:y=0e_1: y=0 上，

ϕ~1(x,0)=(x2−1).\tilde{\boldsymbol\phi}_1(x,0)= \begin{pmatrix} \dfrac{x}{2}\\ -1 \end{pmatrix}.

用外法向 n1=(0,−1)n_1=(0,-1)，可得

ϕ~1⋅n1=1.\tilde{\boldsymbol\phi}_1\cdot n_1 = 1.

注意：底边长度现在是 2，因此总通量为

∫e1ϕ~1⋅n1 ds=1×2=2.\int_{e_1}\tilde{\boldsymbol\phi}_1\cdot n_1\,ds =1\times 2=2.

而参考元对应的自由度是 1。也就是说，直接拉回以后，自由度被放大了一倍。

所以它根本不能作为正确的 H(div)H(\mathrm{div}) 基函数。

11. 散度也会出问题

参考元上

∇x^⋅ϕ^1=∂ξ∂ξ+∂(η−1)∂η=2.\nabla_{\hat x}\cdot \hat{\boldsymbol\phi}_1 = \frac{\partial \xi}{\partial \xi}+ \frac{\partial (\eta-1)}{\partial \eta}=2.

Piola 变换后，理论上应满足

∇x⋅ϕ1=1det⁡J⋅2=1.\nabla_x\cdot \boldsymbol\phi_1 = \frac{1}{\det J}\cdot 2 =1.

而我们确实有

ϕ1=12(x−1y−1),\boldsymbol\phi_1=\frac12\binom{x-1}{y-1},

故

∇x⋅ϕ1=12+12=1.\nabla_x\cdot \boldsymbol\phi_1 = \frac12+\frac12=1.

完全正确。

但如果直接拉回，

ϕ~1=(12x−12yy−1),\tilde{\boldsymbol\phi}_1= \binom{\frac12 x-\frac12 y}{y-1},

则

∇x⋅ϕ~1=12+1=32,\nabla_x\cdot \tilde{\boldsymbol\phi}_1 = \frac12+1=\frac32,

这已经不是正确的变换结果了。

12. 结论

因此，H(div)H(\mathrm{div}) 有限元必须使用 Piola 变换，而不能直接把参考元基函数做普通拉回。原因是：

  1. Piola 变换保持边法向通量自由度不变；
  2. Piola 变换保持散度的正确变换规律；
  3. Piola 变换把参考元上的 RT 空间映成物理元上的 RT 空间；
  4. 只有这样，单元间法向连续性才有意义，H(div)H(\mathrm{div}) 元的定义才成立。

而普通拉回会导致：

  • 法向通量自由度错误；
  • 散度错误；
  • 映射后的函数不再属于物理元上的 RT0 空间。

所以对于 H(div)H(\mathrm{div}) 元，Piola 变换不是一种“方便的选法”，而是由有限元结构本身逼出来的唯一正确选择。

13. 一句话总结

H(div)H(\mathrm{div}) 元要保持的是“法向通量”，不是“向量值本身”；Piola 变换正是保持法向通量和散度结构的变换，因此 H(div)H(\mathrm{div}) 基函数必须通过 Piola 变换定义。

版权归属：Guisong Wu

许可证：署名 4.0 国际 (CC-BY-4.0)