MFEM项目中实现双法向跳跃双线性形式积分器

概述

在有限元分析中，处理间断问题的数值方法经常需要实现特定的面积分器。本文将详细介绍如何在MFEM项目中实现一个用于计算双法向跳跃项的双线性形式积分器，该积分器形式为：

γ₀∫F (1/hF)[u·n][v·n]

其中γ₀是常数，hF表示面的大小度量，u和v属于相同的有限元空间（向量维度vdim=2或3）。

技术实现细节

积分器类设计

我们设计了一个名为NormalJumpIntegrator的类，继承自MFEM的BilinearFormIntegrator基类。该类主要实现面矩阵的组装功能：

class NormalJumpIntegrator : public BilinearFormIntegrator {
public:
   NormalJumpIntegrator(double gamma_0_ = 1.0) : gamma_0(gamma_0_) { }
   
   virtual void AssembleFaceMatrix(const FiniteElement &el1,
                                 const FiniteElement &el2,
                                 FaceElementTransformations &Trans,
                                 DenseMatrix &elmat);
protected:
   double gamma_0;
   Vector shape1, shape2;
   Vector nor;
   
   void AssembleBlock(...);  // 辅助组装函数
};

关键实现要点

1. 法向量处理：通过CalcOrtho函数计算面变换的Jacobian矩阵的正交向量，得到面的法向量nor。

2. 权重计算：考虑了两种情形：

   • 内部面（两个相邻单元）
   • 边界面（仅一个单元）

3. 块矩阵组装：使用AssembleBlock辅助函数组装矩阵的各个块，处理向量值函数的各个分量。

4. 跳跃项处理：通过取负shape函数值来处理法向量方向相反的情况。

实现中的挑战与解决方案

1. 法向量传递问题：最初实现中发现法向量在传递过程中值发生变化，通过移除static关键字和避免传递nor参数解决了此问题。

2. 权重计算：正确计算了从参考面到物理面的变换权重，包括面大小hF的倒数。

3. 矩阵块结构：由于涉及法向量，每个"象限"的矩阵块都是满的，这与某些特殊情形下的稀疏结构不同。

验证与应用

该积分器特别适用于间断有限元方法(DG)中的稳定性项。验证时需要注意：

1. 当应用于连续有限元空间（如H¹）时，理论上应得到零矩阵，这可以作为正确性的初步验证。

2. 对于实际应用，需要根据具体问题调整惩罚参数γ₀的值。

结论

本文详细介绍了在MFEM框架中实现双法向跳跃积分器的关键技术要点。正确实现这类积分器对于处理间断问题和保证数值方法的稳定性至关重要。开发者需要注意法向量的处理、权重计算以及矩阵块结构的特殊性，才能确保积分器的正确性和数值稳定性。