Dual-Mortar对偶基函数与Mortar积分

在上一篇 Mortar 算法介绍 中，我们介绍了 Mortar 方法的基本思想。本文继续讨论 Dual-Mortar 中两个关键问题：

1. Dual basis（对偶基函数）为什么这样构造？
2. 三维非匹配界面上的 Mortar 积分究竟怎么做？

符号约定

• \(\Gamma_s,\Gamma_m\)：slave/master surface；
• \(\Omega_e\)：slave element \(e\) 上参与约束的积分域；
• \(P_{ef}\)：slave element \(e\) 与 master element \(f\) 的 overlap patch；
• \(\mathcal D^e\)：dual basis 构造中的局部对角矩阵；
• \(D\)：最终 Mortar 约束中的 slave-side 矩阵。

本文默认

\[\Omega_e=\bigcup_fP_{ef}. \]

对完整 Tie 界面，它通常就是整个 slave element。若只有部分覆盖或 active region 被截断，则需要额外的边界/基函数处理，不能直接套用下面的“全单元对角化”结论。

1. Mortar 约束问题

首先回顾一下我们在 Mortar算法介绍 中讲的Dual-Mortar 算法。它在总势能中增加了界面约束项：

\[\Pi_c = \int_{\Gamma_s} \boldsymbol\lambda^T \left( \boldsymbol u_s-\boldsymbol u_m \right)d\Gamma \]

其中 \(\boldsymbol\lambda\) 的物理意义是界面牵引力。

形函数插值：

\[\boldsymbol\lambda = \sum_i\hat N_i\boldsymbol\lambda_i, \qquad \boldsymbol u_s = \sum_jN_j^s\boldsymbol d_j^s, \qquad \boldsymbol u_m = \sum_kN_k^m\boldsymbol d_k^m. \]

对拉格朗日乘子取变分，带入形函数积分，最终可得到以下约束方程：

\[D_{ij}d_j^s - C_{ik}d_k^m =0 \]

其中：

\[D_{ij} = \int_{\Gamma_s} \hat{N}_iN_j^s\,d\Gamma \\ C_{ik} = \int_{\Gamma_s} \hat{N}_iN_k^m(\eta(\xi))\,d\Gamma \]

将 slave surface 按 element 分解：

\[\Gamma_s^{\mathrm{active}} = \bigcup_e\Omega_e \]

对单个 slave element，求其局部贡献。

因此全局矩阵可以按 element 和 patch 分解：

\[D_{ij} = \sum_e \sum_f D_{ij}^{ef} \\ C_{ik} = \sum_e \sum_f C_{ik}^{ef} \]

其中：

\[D_{ij}^{ef} = \int_{P_{ef}} \hat{N}_i^eN_j^{s,e} \,d\Gamma \\ C_{ik}^{ef} = \int_{P_{ef}} \hat{N}_i^eN_k^{m,f} \,d\Gamma \]

2. Dual Basis 构造

对 slave element \(e\)，定义局部质量矩阵

\[\boxed{M_{ab}^e=\int_{\Omega_e}N_a^{s,e}N_b^{s,e}\,d\Gamma} \]

以及对角矩阵

\[\boxed{\mathcal D_{ab}^e= \delta_{ab}\int_{\Omega_e}N_a^{s,e}\,d\Gamma}. \]

定义

\[\boxed{A^e=\mathcal D^e(M^e)^{-1}} \]

以及

\[\boxed{\hat N_a^e=\sum_bA_{ab}^eN_b^{s,e}}. \]

于是

\[\begin{aligned} \int_{\Omega_e}\hat N_a^eN_b^{s,e}\,d\Gamma &=\int_{\Omega_e}(A_{ac}^eN_c^{s,e})N_b^{s,e}\,d\Gamma\\ &=A_{ac}^eM_{cb}^e\\ &=(A^eM^e)_{ab}\\ &=\mathcal D_{ab}^e. \end{aligned} \]

因此

\[\boxed{ \int_{\Omega_e}\hat N_a^eN_b^{s,e}\,d\Gamma = \delta_{ab}\int_{\Omega_e}N_a^{s,e}\,d\Gamma } \]

这就是 element-wise biorthogonality（单元级双正交性）。

这里你可能会问：为什么这么构造对偶基函数？我把解释放附录里了。

3. Mortar 积分

对于非匹配三维界面：

• slave/master 节点通常不重合；
• 两个 elements 可能只部分重叠；
• 一个 slave element 可能覆盖多个 master elements；
• slave/master surface 可能不共面；
• 四边形和高阶 surface element 本身甚至可能不是平面。

因此不能简单地拿 slave element 原有的 Gauss points 去计算

\[\int_{\Gamma_s}\hat N_iN_j^m\,d\Gamma. \]

积分区域必须按照 slave/master 的真实 overlap 重新划分。

下面给出一种典型的 3D Mortar integration 流程，基本思想参考 Puso [1]。

3.1 接触搜索与投影

设当前 slave element 为 \(S_k\)，候选 master element 为 \(M_l\)。首先通过接触搜索得到可能发生 overlap 的 master elements。随后为 slave element 或 slave sub-element 构造局部投影平面，并将 slave/master 几何投影到该平面。在投影平面上分别得到：

• slave polygon \(P_s\)；
• master polygon \(P_m\)。

然后求二维多边形交集：

\[\boxed{P=P_s\cap P_m}. \]

这个 \(P\) 就是当前 slave-master element pair 的有效 integration patch。

投影平面的构造本身是一个很 tricky 的问题。对于线性三角形 surface element，单元本身就是平面，因此比较直接；但对于 bilinear quadrilateral、warped quadrilateral 或二阶 surface element，几何并不一定严格共面，此时 projection plane、segmentation 以及几何误差控制都会复杂很多。为了避免深陷太多细节，这里就不深入探讨投影平面的构造了。

3.2 将交叠多边形三角剖分

将 overlap polygon 分解为若干三角形：

\[P=\bigcup_{r=1}^{n_T}T_r. \]

设三角形 \(T_r\) 的三个顶点为

\[\boldsymbol x_0,\qquad \boldsymbol x_1,\qquad \boldsymbol x_2. \]

参考三角形坐标为 \((u,v)\)，则映射为

\[\boldsymbol x(u,v) = (1-u-v)\boldsymbol x_0 + u\boldsymbol x_1 + v\boldsymbol x_2. \]

其面积 Jacobian 为

\[\boxed{ J_T= \left\| (\boldsymbol x_1-\boldsymbol x_0) \times (\boldsymbol x_2-\boldsymbol x_0) \right\| =2A_T }. \]

因此若参考三角形积分点和权重为

\[(u_q,v_q),\qquad \hat w_q, \]

则物理三角形上的积分权重为

\[\boxed{ w_q^{(T)}=J_T\hat w_q }. \]

如果采用 Radon 2D7 等三角形 quadrature rule，只需要把对应的参考积分点和权重代入即可。

3.3 将积分点映射回 slave/master 单元

这里是 Mortar integration 中最容易产生误解的一步。

三角剖分是在 projection plane 上完成的，因此积分点

\[\boldsymbol x_q \]

首先位于投影平面上。

但是我们最终需要计算的是

\[N_a^s(\boldsymbol\xi_s), \qquad N_b^m(\boldsymbol\xi_m), \]

所以必须分别确定这个 quadrature point 对应的 slave/master parametric coordinates：

\[\boxed{ \boldsymbol\xi_{s,q} = \operatorname{InverseMap}_s(\boldsymbol x_q) } \]

\[\boxed{ \boldsymbol\xi_{m,q} = \operatorname{InverseMap}_m(\boldsymbol x_q) }. \]

严格来说，对于曲面或非共面单元，这里的 InverseMap 往往不仅仅是普通的“物理坐标反求等参坐标”，还与前面采用的投影方向和 projection operator 有关。

得到局部坐标以后计算

\[N_{s,a}^{(q)} = N_{s,a}(\boldsymbol\xi_{s,q}), \]

\[N_{m,b}^{(q)} = N_{m,b}(\boldsymbol\xi_{m,q}). \]

因此一个 patch quadrature point 实际同时拥有三套坐标：

projection-plane coordinate: x_q
slave parametric coordinate: ξ_s,q
master parametric coordinate: ξ_m,q

这正是非匹配 Mortar integration 与普通 FEM element integration 的重要区别。

3.4 Slave element 上的 Mortar dual basis

对于 slave element \(e\)，遍历属于它的所有 overlap patches。数值积分得到

\[M_{ac}^e \approx \sum_f \sum_{T\subset P_{ef}} \sum_q w_q^{(T)} N_{s,a}^{(q)}N_{s,c}^{(q)}, \]

以及

\[d_a^e \approx \sum_f \sum_{T\subset P_{ef}} \sum_q w_q^{(T)} N_{s,a}^{(q)}. \]

构造

\[\mathcal D^e = \operatorname{diag}(d_1^e,\ldots,d_n^e), \]

然后

\[\boxed{ A^e=\mathcal D^e(M^e)^{-1} }. \]

因此

\[\boxed{ \hat N_{s,a}^{(q)} = \sum_cA_{ac}^eN_{s,c}^{(q)} }. \]

注意 \(A^e\) 对当前 slave element 的所有 patches 都相同。

3.5 计算 Master-Side Mortar Coefficients

有了 dual basis 后，对每一个 patch 计算

\[G_{aB}^{ef} = \int_{P_{ef}} \hat N_{s,a}^eN_{m,B}^f \,dS. \]

数值积分为

\[\boxed{ G_{aB}^{ef} \approx \sum_{T\subset P_{ef}} \sum_q w_q^{(T)} \hat N_{s,a}^e(\boldsymbol\xi_{s,q}) N_{m,B}^f(\boldsymbol\xi_{m,q}) }. \]

然后对与 slave element \(e\) 相交的所有 master elements 求和：

\[G_{aB}^e = \sum_fG_{aB}^{ef}. \]

这里 \(B\) 是 master surface 的全局节点编号，因此实现时需要把每个 master element 的局部列贡献 scatter 到全局 master DOFs。

3.6 装配到 Slave Node

对于 slave surface 上的全局节点 \(A\)，它可能属于多个 slave elements。

局部对角系数装配为

\[\boxed{ d_A = \sum_{e\ni A}d_A^e }. \]

master coupling coefficients 装配为

\[\boxed{ G_{AB} = \sum_{e\ni A}G_{AB}^e }. \]

因此全局约束的第 \(A\) 行为

\[\boxed{ d_A\boldsymbol u_{s,A} - \sum_BG_{AB}\boldsymbol u_{m,B} = \boldsymbol0 }. \]

由于 \(d_A\) 是标量，可以直接得到

\[\boxed{ \boldsymbol u_{s,A} = \sum_B \frac{G_{AB}}{d_A} \boldsymbol u_{m,B} }. \]

展开就是

\[\boxed{ \boldsymbol u_{s,A} = \sum_B \frac{ \displaystyle \sum_{e\ni A} \int_{\Omega_e} \hat N_{s,A}^eN_{m,B}\,dS }{ \displaystyle \sum_{e\ni A} \int_{\Omega_e} N_{s,A}^e\,dS } \boldsymbol u_{m,B} }. \]

这就是 Dual-Mortar 最终非常漂亮的结果：

slave node 的位移可以写成一组 master nodal displacements 的加权组合，而不需要求解一个全局的 \(D^{-1}\)。

3.7 总流程

把前面的内容串起来，算法可以概括为：

for each slave element e:

    1. 搜索 candidate master elements

    2. 对每个 candidate master element f:
           构造 projection plane
           投影 slave/master element
           polygon clipping
           得到 overlap patch P_ef
           triangulate P_ef
           生成 quadrature points

    3. 利用当前 slave element 的所有 patches:
           accumulate M_e
           accumulate d_e

    4. 构造:
           Dcal_e = diag(d_e)
           A_e = Dcal_e * inverse(M_e)

    5. 再遍历 patches:
           Nhat_s = A_e * N_s
           G_e += ∫ Nhat_s^T N_m dS

    6. assemble:
           d_global += d_e
           G_global += G_e

最终：

for each slave node A:
    u_s[A] = Σ_B (G[A,B] / d[A]) * u_m[B]

从矩阵角度看，就是

\[D\boldsymbol u_s-C\boldsymbol u_m=0, \]

而 Dual-Mortar 通过 biorthogonal multiplier basis 使

\[D=\operatorname{diag}(d_1,d_2,\ldots). \]

因此

\[\boxed{ \boldsymbol u_s = D^{-1}C\boldsymbol u_m } \]

中的 \(D^{-1}\) 只是逐行除以对角元素。

4. Consistency Check

如果完整覆盖、几何映射和数值积分都正确，那么

\[\boxed{ \sum_Bw_{AB}=1 }. \]

因此它可以作为 Mortar implementation 的一个非常实用的单元测试 / consistency check。

5. 总结

Dual-Mortar 的核心其实可以浓缩成两件事。

第一件事是 代数上的对角化：通过施加 element-wise biorthogonality 使得全局 slave-side Mortar matrix \(D\) 成为对角矩阵。

第二件事是 几何上的非匹配积分：

\[\boxed{ \text{search} \rightarrow \text{projection} \rightarrow \text{clipping} \rightarrow \text{triangulation} \rightarrow \text{quadrature} \rightarrow \text{inverse mapping} \rightarrow \text{assembly} } \]

前者解决“如何让约束容易消元”，后者解决“两个不匹配网格之间究竟在哪里、怎么积分”。

参考文献

[1] Puso, Michael A. ‘A 3D Mortar Method for Solid Mechanics’. International Journal for Numerical Methods in Engineering 59, no. 3 (2004): 315–36. https://doi.org/10.1002/nme.865.

附录

对偶基函数

我们希望寻找 \(\hat N_i\)，满足：

\[\boxed{ \langle\hat N_i,N_j\rangle = \int_{\Omega_e} \hat N_iN_j\,d\Gamma = c_i\delta_{ij} } \]

这里先不要管 \(c_i\) 为什么这么选。它表达的核心意思就是：\(i\neq j \quad\Rightarrow\quad \int\hat N_iN_j\,d\Gamma=0.\)

也就是：

\(\hat N_i\) 在积分意义下只和对应的 \(N_i\) 耦合，与其他 \(N_j\) 不耦合。

这就是双正交性 biorthogonality。

那么 \(\hat N_i\) 从哪里找？我们引入个限制，要求 dual shape functions 位于 slave shape functions 张成的空间中：

\[\boxed{ \hat N_i = \sum_a A_{ia}N_a } \]

为了最后不计算一个大型稠密矩阵的逆，我们主动施加更强的 element-wise biorthogonality，即每一个 element 都要求：

\[\int_{\Omega_e} \hat N_a^eN_b^e\,d\Gamma = c_a^e\delta_{ab} \]

然后利用有限元装配 \(\int_{\Gamma_s} = \sum_e\int_{\Omega_e}\) ，这样得到的全局矩阵仍然是对角的。

我们举一个简单的例子解释一下上面的公式。假设 slave surface 是：

             Γs

      e1              e2
1 ----------- 2 ----------- 3

单元 \(e_1\) 有局部节点：\((1,2)\)，单元 \(e_2\)：\((2,3).\)

我们分别构造：\(D^{e_1} = \begin{bmatrix} d_1^{e_1}&0\\ 0&d_2^{e_1} \end{bmatrix}\) 和 \(D^{e_2} = \begin{bmatrix} d_2^{e_2}&0\\ 0&d_3^{e_2} \end{bmatrix}.\)

然后像普通 FEM matrix 一样 assemble。最终：

\[D= \begin{bmatrix} d_1^{e_1}&0&0\\ 0&d_2^{e_1}+d_2^{e_2}&0\\ 0&0&d_3^{e_2} \end{bmatrix} \]

这里有一个微妙的问题：如果 \(\hat N_2^{e_1}\) 是在 \(e_1\) 上单独构造的，\(\hat N_2^{e_2}\) 又是在 \(e_2\) 上单独构造的，那么共享节点 2 的 dual shape function 到底是哪一个？

这就是 element-wise dual basis 的一个重要特点。它可以理解为一个分片定义（piecewise-defined）的 dual function：

\[\hat N_2(\mathbf x) = \begin{cases} \hat N_2^{e_1}(\mathbf x), & \mathbf x\in\Omega_{e_1}, \\[4pt] \hat N_2^{e_2}(\mathbf x), & \mathbf x\in\Omega_{e_2}, \\ 0,&\text{otherwise}. \end{cases} \]

它不一定像普通 \(C^0\) FEM shape function 那样，在 element boundary 上具有相同的连续性要求。这是允许的，因为：

Lagrange multiplier / dual space 不一定需要和 displacement space 具有相同的连续性。

它主要承担的是 weak constraint 的测试/乘子空间角色。