第26章 积分点、雅可比矩阵、基函数与内部边/面

---

26.1 学习目标

• 理解有限元核心计算数据流：从参考单元积分点（Integration Points）到物理空间雅可比矩阵（Jacobian Matrix）再到基函数（Basis Functions）求值
• 掌握 mesh::getIntegrationPoints() 获取高斯积分规则（Gauss-Legendre）的局部坐标与权重
• 掌握 mesh::getJacobians() 批量计算所有单元的雅可比矩阵、行列式及物理坐标——理解按列返回的内存布局
• 掌握 mesh::getBasisFunctions() 获取 Lagrange 基函数值及其梯度（GradLagrange），以及高阶分层基函数（H1Legendre / HcurlLegendre）
• 掌握 mesh::getElementEdgeNodes() / getElementFaceNodes() — 提取每个单元边/面的节点列表
• 掌握 mesh::createEdges() / createFaces() 创建全局唯一边/面，及 getAllEdges() / getAllFaces() 一次性获取
• 学会在离散实体上创建基于内部面的新单元——DG（Discontinuous Galerkin）方法前处理的关键步骤
• 了解 mesh::getElementProperties() 与 mesh::getIntegrationData() 等辅助 API 获取单元元信息与综合积分数据

---

26.2 核心概念说明

26.2.1 从参考单元到物理单元的映射

在标准有限元方法中，刚度矩阵和质量矩阵的组装遵循以下流程：

参考单元 (Reference Element)                      物理单元 (Physical Element)
┌─────────────────────────┐                       ┌──────────────────────┐
│  积分点: (u, v, w)       │                       │  积分点: (x, y, z)    │
│  权重:   w_i            │                       │                       │
│  基函数: N_j(u, v, w)   │  ── 雅可比映射 ──>    │  基函数: N_j(x, y, z) │
│  梯度: ∂N_j/∂(u,v,w)    │                       │  梯度:  ∂N_j/∂(x,y,z) │
└─────────────────────────┘                       └──────────────────────┘

雅可比矩阵 J 是从参考坐标 (u, v, w) 到物理坐标 (x, y, z) 的偏导数矩阵：

        [ ∂x/∂u  ∂x/∂v  ∂x/∂w ]
    J = [ ∂y/∂u  ∂y/∂v  ∂y/∂w ]          (3x3 矩阵，按列存储)
        [ ∂z/∂u  ∂z/∂v  ∂z/∂w ]

物理空间梯度通过链式法则计算：∂N/∂x = J^{-T} * ∂N/∂u，其中 det(J) 用于积分变量变换 dx dy dz = |det(J)| du dv dw。

Gmsh 通过三个核心 API 提供这些数据： 1. getIntegrationPoints() — 参考单元中的积分点位置和权重 2. getJacobians() — 各积分点处的雅可比矩阵、行列式和物理坐标 3. getBasisFunctions() — 各积分点处的基函数值和梯度

26.2.2 积分规则家族（Integration Rule Families）

Gmsh 支持多种积分规则，通过 integrationType 字符串参数指定：

规则家族	integrationType 示例	说明
Gauss	"Gauss2", "Gauss4"	经济型高斯积分（点数最少），不可用时回退到 CompositeGauss
CompositeGauss	"CompositeGauss4"	张量积规则（基于一维 Gauss-Legendre），始终可用
GaussLegendre	"GaussLegendre4"	标准 Gauss-Legendre 积分（一维）

选择依据： - 对四边形/六面体，"Gauss" 使用点最少的高斯规则 - 对三角形/四面体，"Gauss" 使用对称经济型规则 - 如需确定积分点都在参考单元内部，推荐 "CompositeGauss"（"Gauss" 高阶时个别点可能落在参考单元外）

26.2.3 基函数空间类型（Function Space Types）

getBasisFunctions() 支持多种函数空间，通过 functionSpaceType 参数指定：

函数空间	示例参数	含义	每个积分点返回的分量数
Lagrange	"Lagrange"	等参 Lagrange 基函数值	1（标量值）
LagrangeN	"Lagrange2"	N 阶 Lagrange 基函数值	1
GradLagrange	"GradLagrange"	等参 Lagrange 基函数梯度	3（∂/∂u, ∂/∂v, ∂/∂w）
GradLagrangeN	"GradLagrange2"	N 阶 Lagrange 基函数梯度	3
H1LegendreN	"H1Legendre3"	N 阶分层 H1 Legendre 基函数	1
GradH1LegendreN	"GradH1Legendre3"	N 阶分层 H1 Legendre 梯度	3
HcurlLegendreN	"HcurlLegendre3"	N 阶分层 H(curl) Legendre 基函数	3（矢量）
CurlHcurlLegendreN	"CurlHcurlLegendre3"	N 阶分层 H(curl) Legendre 旋度	3

分层基函数（Hierarchical Basis Functions）允许 p-细化（p-refinement）而不改变低阶基函数，对 hp-自适应方法特别有利。

26.2.4 内部边/面与 DG 方法

在 DG（Discontinuous Galerkin）方法中，单元间通过内部面（Internal Face）的数值通量（Numerical Flux）进行耦合。Gmsh 的网格边/面体系专门为此类需求设计：

• 网格边（Mesh Edge）：一阶线单元的两个节点唯一确定一条边，Gmsh 可为其分配全局唯一边标签，并记录方向（Orientation）
• 网格面（Mesh Face）：三角形面（3 个节点）或四边形面（4 个节点），同样可分配全局唯一标签
• 内部边/面：指被两个或多个单元共享的边/面——区别于位于模型边界上的边界边/面

Gmsh 工作流： 1. createEdges() / createFaces() — 创建全局唯一边/面 2. getElementEdgeNodes() / getElementFaceNodes() — 获取每个单元的边/面节点 3. getEdges() / getFaces() — 通过节点查询边/面标签 4. 在离散实体上新建单元——将内部面绑定到新的低维单元上，便于积分和插值

---

26.3 C++ 代码逐段讲解

26.3.1 头文件、初始化与辅助函数

#include <iostream>
#include <set>
#include <map>
#include <gmsh.h>

int main(int argc, char **argv)
{
  gmsh::initialize(argc, argv);
  gmsh::model::add("ch26");

  // 辅助打印函数：输出向量内容，mult 控制每几个数一组的标签
  auto pp = [](const std::string &label,
              const std::vector<double> &v, int mult)
  {
    std::cout << " * " << v.size() / mult << " " << label << ": ";
    for(auto c : v) std::cout << c << " ";
    std::cout << "\n";
  };

解释：pp 是一个 lambda 辅助函数，用于格式化输出向量内容。mult 参数指示每组有几个分量——例如 localCoords 每组 3 个分量 (u, v, w)，则 mult = 3；Jacobian determinants 每组 1 个值，则 mult = 1。

---

26.3.2 创建几何与网格

  // 创建 2D 矩形域：宽 1.0，高 0.1
  gmsh::model::occ::addRectangle(0, 0, 0, 1, 0.1);
  gmsh::model::occ::synchronize();

  // 使用 transfinite 算法生成结构化四边形网格
  gmsh::model::mesh::setTransfiniteAutomatic();

  // 设置一阶单元 (elementOrder = 1) 并生成 2D 网格
  int elementOrder = 1, interpolationOrder = 2;
  gmsh::model::mesh::setOrder(elementOrder);
  gmsh::model::mesh::generate(2);

解释：创建一个 1.0 x 0.1 的矩形，用 transfinite 自动算法生成结构化四边形网格。elementOrder = 1 表示一阶单元（4 节点四边形）；interpolationOrder = 2 将用于选择能精确积分二次多项式的积分规则。

---

26.3.3 遍历单元类型并获取单元属性

  // 遍历网格中所有单元类型
  std::vector<int> elementTypes;
  gmsh::model::mesh::getElementTypes(elementTypes);

  for(auto t : elementTypes) {
    // 获取该单元类型的属性
    std::string elementName;
    int dim, order, numNodes, numPrimNodes;
    std::vector<double> localNodeCoord;
    gmsh::model::mesh::getElementProperties(
        t, elementName, dim, order, numNodes,
        localNodeCoord, numPrimNodes);
    std::cout << "\n** " << elementName << " **\n\n";

解释：getElementTypes() 返回当前网格中所有 MSH 单元类型的整数标识（如 3 = 4 节点四边形 Quadrangle 4）。getElementProperties() 返回该类型的： - elementName：名称字符串（如 "Quadrangle 4"） - dim：几何维度 - order：几何阶数（对一阶单元 order = 1） - numNodes：单元总节点数 - localNodeCoord：参考单元中所有节点的局部坐标，长度 = dim * numNodes - numPrimNodes：（primary nodes）一阶顶点节点数

---

26.3.4 获取积分点（Integration Points）

    // -- A. 积分点 --
    // 获取能精确积分 "interpolationOrder" 阶多项式的 Gauss 积分规则
    std::vector<double> localCoords, weights;
    gmsh::model::mesh::getIntegrationPoints(
        t, "Gauss" + std::to_string(interpolationOrder),
        localCoords, weights);
    pp("integration points to integrate order " +
       std::to_string(interpolationOrder) +
       " polynomials", localCoords, 3);

解释：getIntegrationPoints(elementType, integrationType, localCoord, weights) 返回：

• localCoord：所有积分点在参考单元中的局部坐标，按 [g1u, g1v, g1w, g2u, g2v, g2w, ...] 顺序拼接。二维单元中 w 坐标恒为 0。
• weights：对应的积分权重 [w1, w2, ..., wG]。

这里使用 "Gauss2" 规则，对于一阶四边形单元（4 节点，bilinear 插值），Gauss 积分需要 2 x 2 = 4 个积分点来精确积分二次多项式（由于张量积规则）。每条坐标轴上 2 个 Gauss 点，局部坐标为 ±1/√3 ≈ ±0.57735。

典型输出示例（Quadrangle 4, "Gauss2"）：

 * 4 integration points to integrate order 2 polynomials:
   -0.57735 -0.57735 0  0.57735 -0.57735 0
   0.57735  0.57735 0 -0.57735  0.57735 0

---

26.3.5 获取基函数值及其梯度

    // -- C. 基函数 --
    // 获取 Lagrange 基函数在积分点上的值
    int numComponents, numOrientations;
    std::vector<double> basisFunctions;

    gmsh::model::mesh::getBasisFunctions(
        t, localCoords, "Lagrange",
        numComponents, basisFunctions, numOrientations);
    pp("basis functions at integration points",
       basisFunctions, 1);

    // 获取 Lagrange 基函数的梯度 (∂N/∂u, ∂N/∂v, ∂N/∂w) 在积分点上的值
    gmsh::model::mesh::getBasisFunctions(
        t, localCoords, "GradLagrange",
        numComponents, basisFunctions, numOrientations);
    pp("basis function gradients at integration points",
       basisFunctions, 3);

解释：getBasisFunctions(elementType, localCoord, functionSpaceType, numComponents, basisFunctions, numOrientations)：

• "Lagrange"：返回每个积分点上每个基函数的标量值。对于具有 N 个节点的单元，每个积分点返回 N 个值。数据按 [g1N1, g1N2, ..., g1NN, g2N1, g2N2, ...] 排列。
• "GradLagrange"：返回每个积分点上每个基函数的梯度向量 (∂N/∂u, ∂N/∂v, ∂N/∂w)。每个基函数每个积分点贡献 3 个值。数据按 [g1N1_du, g1N1_dv, g1N1_dw, g1N2_du, ..., gGN_NN_dw] 排列。
• numComponents：每个基函数在每个积分点返回的分量数（标量=1，梯度=3）
• numOrientations：基函数方向总数（Lagrange 基函数为 1——各向同性）

对于一阶四边形单元在 4 个 Gauss 积分点上： - "Lagrange" 返回 4x4 = 16 个值 - "GradLagrange" 返回 4x4x3 = 48 个值

高阶基函数示例：若调用 "Lagrange2" 则返回二阶 Lagrange 基函数值，适用于 p=2 的有限元空间；"H1Legendre3" 则返回三阶分层 H1 Legendre 基函数，适用于 hp-自适应方法。

---

26.3.6 获取雅可比矩阵

    // -- B. 雅可比矩阵 --
    // 计算所有该类型单元的雅可比矩阵、行列式及物理坐标
    std::vector<double> jacobians, determinants, coords;
    gmsh::model::mesh::getJacobians(
        t, localCoords, jacobians, determinants, coords);
    pp("Jacobian determinants at integration points",
       determinants, 1);

解释：getJacobians(elementType, localCoord, jacobians, determinants, coord, tag, task, numTasks) 批量计算所有属于该类型的单元在各个积分点处的雅可比数据：

• jacobians：3x3 雅可比矩阵按列存储。每个单元每个积分点存储 9 个值：[Jxu, Jyu, Jzu, Jxv, Jyv, Jzv, Jxw, Jyw, Jzw]，其中 Jxu = ∂x/∂u, Jyu = ∂y/∂u 等。数据排列方式为： [e1_g1_col1, e1_g1_col2, e1_g1_col3, e1_g2_col1, ..., e1_gG_col3, e2_g1_col1, ...]
• determinants：每个单元每个积分点的 det(J) 值，[e1g1, e1g2, ..., e1gG, e2g1, ...]
• coord：每个积分点在物理空间中的 (x, y, z) 坐标，[e1g1_x, e1g1_y, e1g1_z, e1g2_x, ...]

按列存储说明：jacobians 中存储的 9 个值对应矩阵的各列（而非各行）：

  offset 0: J(0,0)=∂x/∂u    offset 3: J(0,1)=∂x/∂v    offset 6: J(0,2)=∂x/∂w
  offset 1: J(1,0)=∂y/∂u    offset 4: J(1,1)=∂y/∂v    offset 7: J(1,2)=∂y/∂w
  offset 2: J(2,0)=∂z/∂u    offset 5: J(2,1)=∂z/∂v    offset 8: J(2,2)=∂z/∂w

这是 Gmsh 的内部存储惯例，与许多有限元库的按行存储不同，在组装刚度矩阵时需注意索引换算。

并行计算支持：task 和 numTasks 参数用于多线程并行计算（如 OpenMP），将总工作量拆分为 numTasks 份，每次调用仅计算 task 编号对应的那一份。使用前需通过 preallocateJacobians() 预分配输出向量。

单元素版本：mesh::getJacobian(elementTag, localCoord, jacobians, determinants, coord) 获取单个单元（通过标签指定）的雅可比数据。

---

26.3.7 第二部分：创建体网格并获取边/面数据（来自 x7.cpp）

  // ===== 第二部分：内部边/面 =====
  // 创建 3D 模型并生成四面体网格
  gmsh::model::add("ch26_3d");
  gmsh::model::occ::addBox(0, 0, 0, 1, 1, 1);
  gmsh::model::occ::synchronize();
  gmsh::option::setNumber("Mesh.MeshSizeMin", 0.2);
  gmsh::model::mesh::generate(3);

解释：创建一个单位立方体，设置最小网格尺寸为 0.2（粗网格便于观察），生成三维四面体网格。这是演示内部边/面操作的标准场景——体网格内部面是两个相邻四面体的共享三角形面。

---

26.3.8 获取每个单元的边节点和面节点

  // 获取一阶四面体的 MSH 类型标识
  int tetType = gmsh::model::mesh::getElementType("tetrahedron", 1);

  // 获取所有四面体单元的边节点列表
  std::vector<std::size_t> edgeNodes, faceNodes;
  gmsh::model::mesh::getElementEdgeNodes(tetType, edgeNodes);

  // 获取所有四面体单元的面节点列表（三角形面，faceType = 3）
  gmsh::model::mesh::getElementFaceNodes(tetType, 3, faceNodes);

解释： - getElementType(familyName, order)：根据单元系列名称和阶数返回 MSH 类型整数标识。"tetrahedron" 对应 4 节点四面体（MSH 类型 4）。 - getElementEdgeNodes(elementType, nodeTags)：返回所有该类型单元的边节点。每个四面体有 6 条边，每条边由 2 个节点（对于一阶单元）或更多（对于高阶单元）定义。数据按单元顺序排列：[e1_edge1_n1, e1_edge1_n2, e1_edge2_n1, ...]。若设置 primary = true，则仅返回边的起点和终点节点（忽略高阶中间节点）。 - getElementFaceNodes(elementType, faceType, nodeTags)：类似地获取面节点。faceType=3 对应三角形面，faceType=4 对应四边形面。每个四面体有 4 个三角形面。

关键点：边节点和面节点的返回顺序与 getElements() / getElementsByType() 中的单元顺序完全一致，因此可以通过位置索引直接将边/面关联到所属单元。

---

26.3.9 创建全局唯一边/面并建立单元关联

  // 创建全局唯一的网格边和面（内部识别，分配标签）
  gmsh::model::mesh::createEdges();
  gmsh::model::mesh::createFaces();

  // 通过节点列表查询边和面的全局标签及方向
  std::vector<std::size_t> edgeTags, faceTags;
  std::vector<int> edgeOrientations, faceOrientations;
  gmsh::model::mesh::getEdges(edgeNodes, edgeTags, edgeOrientations);
  gmsh::model::mesh::getFaces(3, faceNodes, faceTags, faceOrientations);

解释： - createEdges() / createFaces()：扫描网格，为每条唯一边和每个唯一面分配一个全局标签。一条边被多个单元共享时，只分配一个标签，从而可以识别内部边/面（标签被多个单元引用即内部，被单个单元引用即边界）。 - getEdges(nodeTags, edgeTags, edgeOrientations)：输入边节点列表（2N 个节点标识，N = 边数），返回每条边的全局标签和方向。方向 1 表示节点顺序与参考正方向一致，-1 表示相反。 - getFaces(faceType, nodeTags, faceTags, faceOrientations)：同上，作用于面。faceType 必须与面节点数匹配（3 = 三角形，4 = 四边形）。

---

26.3.10 建立边/面到单元的映射

  // 获取所有四面体单元标签和节点，建立边/面到单元的映射
  std::vector<std::size_t> elementTags, elementNodeTags;
  gmsh::model::mesh::getElementsByType(
      tetType, elementTags, elementNodeTags);

  // 边到单元的映射（每个四面体 6 条边）
  std::map<std::size_t, std::vector<std::size_t>> edges2elems;
  for(std::size_t i = 0; i < edgeTags.size(); i++)
    edges2elems[edgeTags[i]].push_back(elementTags[i / 6]);

  // 面到单元的映射（每个四面体 4 个面）
  std::map<std::size_t, std::vector<std::size_t>> faces2elems;
  for(std::size_t i = 0; i < faceTags.size(); i++)
    faces2elems[faceTags[i]].push_back(elementTags[i / 4]);

  // 打印内部面信息
  std::cout << "\n** Internal faces (shared by 2 tetrahedra) **\n";
  for(auto &kv : faces2elems) {
    if(kv.second.size() == 2) {
      std::cout << "Face " << kv.first
                << " connects tet " << kv.second[0]
                << " and tet " << kv.second[1] << "\n";
    }
  }

解释：这是 DG 方法中构建单元邻接关系的关键步骤： 1. getElementsByType() 按单元类型获取单元标签和节点，返回顺序与边/面节点相同 2. 整数除法 i / 6（或 i / 4）将边/面索引映射回单元索引——因为每个四面体贡献固定数量的边和面 3. edges2elems 和 faces2elems 映射记录了每条边/面被哪些单元共享 4. 共享计数为 2 的边/面是内部边/面（Internal Edges/Faces），为 1 的是边界边/面

---

26.3.11 在离散实体上创建基于内部面的新单元

  // 创建离散面实体，将共享面绑定到新单元（DG 方法关键步骤）
  int s = gmsh::model::addDiscreteEntity(2);

  // 收集唯一面标签，为每个唯一内部面创建二维三角形单元
  std::set<std::size_t> uniqueFaceTags;
  std::vector<std::size_t> triTags, triNodes;
  std::size_t maxTag;
  gmsh::model::mesh::getMaxElementTag(maxTag);

  for(std::size_t i = 0; i < faceTags.size(); i++) {
    if(uniqueFaceTags.find(faceTags[i]) == uniqueFaceTags.end()) {
      uniqueFaceTags.insert(faceTags[i]);
      // 基于面标签创建新的单元标签（避免与已有标签冲突）
      triTags.push_back(faceTags[i] + maxTag);
      // 提取该面的三个节点
      triNodes.push_back(faceNodes[3 * i]);
      triNodes.push_back(faceNodes[3 * i + 1]);
      triNodes.push_back(faceNodes[3 * i + 2]);
    }
  }

  int triType = gmsh::model::mesh::getElementType("triangle", 1);
  gmsh::model::mesh::addElementsByType(
      s, triType, triTags, triNodes);

解释：这段代码是 DG 前处理的核心操作： 1. addDiscreteEntity(2) 在模型中添加一个二维离散实体——它没有参数化几何，仅由网格数据定义 2. getMaxElementTag() 获取当前最大单元标签，用于为新单元分配不冲突的标签 3. 遍历所有面，跳过重复的标签（uniqueFaceTags 集合去重），为每个唯一面创建一个二维三角形单元 4. addElementsByType(entityTag, elementType, elementTags, nodeTags) 将新单元添加到离散实体上

为什么这样做：在 DG 方法中，内部面上的数值积分（计算数值通量）是核心操作。将内部面创建为独立的低维单元后，便可以利用 26.3.4-26.3.6 节的 API（getIntegrationPoints()、getJacobians()、getBasisFunctions()）直接在面上进行积分，如同处理普通单元一样。

---

26.3.12 邻接信息输出与一次性 API

  // 打印三角形单元与其相邻四面体的关系
  for(auto tTag : triTags) {
    std::cout << "triangle " << tTag
              << " connected to tetrahedra:";
    for(auto tt : faces2elems[tTag - maxTag])
      std::cout << " " << tt;
    std::cout << "\n";
  }

  // 一次性获取所有边/面的便捷 API
  std::vector<std::size_t> allEdgeTags, allEdgeNodes;
  gmsh::model::mesh::getAllEdges(allEdgeTags, allEdgeNodes);

  std::vector<std::size_t> allFaceTags, allFaceNodes;
  gmsh::model::mesh::getAllFaces(3, allFaceTags, allFaceNodes);

  std::cout << "Total unique edges: " << allEdgeTags.size() << "\n";
  std::cout << "Total unique faces: " << allFaceTags.size() << "\n";

解释： - 通过 tTag - maxTag 反查原始面标签，即可从 faces2elems 映射中获得该面所属的四面体单元——这在 DG 中用于确定通量计算的左/右单元 - getAllEdges(edgeTags, edgeNodes) 和 getAllFaces(faceType, faceTags, faceNodes) 是一步到位的便捷 API，直接返回网格中所有唯一边/面的标签和节点，无需事先调用 getElementEdgeNodes() 或手动收集。适用于仅需要边/面拓扑信息而不关心单元关联的场景

---

26.3.13 GUI 与清理

  // 启动 GUI（可通过 -nopopup 跳过）
  std::set<std::string> args(argv, argv + argc);
  if(!args.count("-nopopup")) gmsh::fltk::run();

  gmsh::finalize();
  return 0;
}

解释：gmsh::fltk::run() 启动交互式 GUI 以可视化结果。支持 -nopopup 命令行参数跳过 GUI（适用于批处理/自动化脚本）。

---

26.4 完整可运行代码

// =============================================================================
//  Chapter 26: 积分点、雅可比矩阵、基函数与内部边/面
//
//  功能: 演示 Gmsh 有限元核心数据 API：
//        (1) 积分点坐标与权重
//        (2) 雅可比矩阵、行列式及物理坐标
//        (3) Lagrange 基函数值及梯度
//        (4) 网格内部边/面的创建、查询与 DG 前处理
//
//  编译 (需链接 Gmsh SDK):
//    g++ -o ch26_fe_data ch26_fe_data.cpp -lgmsh -std=c++11
//
//  运行:
//    ./ch26_fe_data              (带 GUI 查看)
//    ./ch26_fe_data -nopopup     (仅文本输出)
// =============================================================================

#include <iostream>
#include <set>
#include <map>
#include <gmsh.h>

int main(int argc, char **argv)
{
  gmsh::initialize(argc, argv);

  // 辅助函数：格式化打印向量内容
  auto pp = [](const std::string &label,
              const std::vector<double> &v, int mult)
  {
    std::cout << " * " << v.size() / mult << " " << label << ": ";
    for(auto c : v) std::cout << c << " ";
    std::cout << "\n\n";
  };

  // ===== Part 1: 2D 积分点、雅可比、基函数 =====
  {
    gmsh::model::add("ch26_2d");
    gmsh::model::occ::addRectangle(0, 0, 0, 1, 0.1);
    gmsh::model::occ::synchronize();
    gmsh::model::mesh::setTransfiniteAutomatic();

    int elementOrder = 1, interpolationOrder = 2;
    gmsh::model::mesh::setOrder(elementOrder);
    gmsh::model::mesh::generate(2);

    std::vector<int> elementTypes;
    gmsh::model::mesh::getElementTypes(elementTypes);

    for(auto t : elementTypes) {
      std::string elemName;
      int dim, order, numNodes, numPrimNodes;
      std::vector<double> localNodeCoord;
      gmsh::model::mesh::getElementProperties(
          t, elemName, dim, order, numNodes,
          localNodeCoord, numPrimNodes);
      std::cout << "\n========== " << elemName
                << " (type " << t << ") ==========\n\n";

      // ---- A. 积分点 ----
      std::vector<double> localCoords, weights;
      gmsh::model::mesh::getIntegrationPoints(
          t, "Gauss" + std::to_string(interpolationOrder),
          localCoords, weights);
      pp("integration points (u,v,w per point)",
         localCoords, 3);
      std::cout << "     weights: ";
      for(auto w : weights) std::cout << w << " ";
      std::cout << "\n\n";

      // ---- C. 基函数值 ----
      int numComp, numOri;
      std::vector<double> basisFunc;
      gmsh::model::mesh::getBasisFunctions(
          t, localCoords, "Lagrange",
          numComp, basisFunc, numOri);
      pp("Lagrange basis function values at integration points",
         basisFunc, 1);

      // ---- C. 基函数梯度 ----
      gmsh::model::mesh::getBasisFunctions(
          t, localCoords, "GradLagrange",
          numComp, basisFunc, numOri);
      pp("GradLagrange (dN/du, dN/dv, dN/dw) at integration points",
         basisFunc, 3);

      // ---- B. 雅可比矩阵 ----
      std::vector<double> jacobians, determinants, coords;
      gmsh::model::mesh::getJacobians(
          t, localCoords, jacobians, determinants, coords);
      pp("Jacobian determinants at integration points",
         determinants, 1);
    }
  }

  // ===== Part 2: 3D 内部边/面 =====
  {
    gmsh::model::add("ch26_3d");
    gmsh::model::occ::addBox(0, 0, 0, 1, 1, 1);
    gmsh::model::occ::synchronize();
    gmsh::option::setNumber("Mesh.MeshSizeMin", 0.2);
    gmsh::model::mesh::generate(3);

    int tetType = gmsh::model::mesh::getElementType("tetrahedron", 1);

    // 获取每个单元的边节点和面节点
    std::vector<std::size_t> edgeNodes, faceNodes;
    gmsh::model::mesh::getElementEdgeNodes(tetType, edgeNodes);
    gmsh::model::mesh::getElementFaceNodes(tetType, 3, faceNodes);

    // 创建全局唯一边/面
    gmsh::model::mesh::createEdges();
    gmsh::model::mesh::createFaces();

    // 查询边/面标签
    std::vector<std::size_t> edgeTags, faceTags;
    std::vector<int> edgeOrient, faceOrient;
    gmsh::model::mesh::getEdges(edgeNodes, edgeTags, edgeOrient);
    gmsh::model::mesh::getFaces(3, faceNodes, faceTags, faceOrient);

    // 获取单元数据，建立映射
    std::vector<std::size_t> elemTags, elemNodeTags;
    gmsh::model::mesh::getElementsByType(
        tetType, elemTags, elemNodeTags);

    std::map<std::size_t, std::vector<std::size_t>> faces2elems;
    for(std::size_t i = 0; i < faceTags.size(); i++)
      faces2elems[faceTags[i]].push_back(elemTags[i / 4]);

    // 打印内部面
    std::cout << "===== Internal Faces =====\n";
    for(auto &kv : faces2elems) {
      if(kv.second.size() == 2) {
        std::cout << "  Face " << kv.first
                  << " : tet " << kv.second[0]
                  << " <-> tet " << kv.second[1] << "\n";
      }
    }

    // 在离散实体上创建基于内部面的新单元
    int s = gmsh::model::addDiscreteEntity(2);
    std::set<std::size_t> uniqueFaceTags;
    std::vector<std::size_t> triTags, triNodes;
    std::size_t maxTag;
    gmsh::model::mesh::getMaxElementTag(maxTag);

    for(std::size_t i = 0; i < faceTags.size(); i++) {
      if(uniqueFaceTags.find(faceTags[i]) == uniqueFaceTags.end()) {
        uniqueFaceTags.insert(faceTags[i]);
        triTags.push_back(faceTags[i] + maxTag);
        triNodes.push_back(faceNodes[3 * i]);
        triNodes.push_back(faceNodes[3 * i + 1]);
        triNodes.push_back(faceNodes[3 * i + 2]);
      }
    }

    int triType = gmsh::model::mesh::getElementType("triangle", 1);
    gmsh::model::mesh::addElementsByType(s, triType, triTags, triNodes);

    // 一次性获取所有边/面
    std::vector<std::size_t> allE, allEN, allF, allFN;
    gmsh::model::mesh::getAllEdges(allE, allEN);
    gmsh::model::mesh::getAllFaces(3, allF, allFN);
    std::cout << "\nUnique edges: " << allE.size()
              << ", Unique faces: " << allF.size() << "\n";
  }

  // GUI
  std::set<std::string> args(argv, argv + argc);
  if(!args.count("-nopopup")) gmsh::fltk::run();

  gmsh::finalize();
  return 0;
}

---

26.5 关键 API 速查表

26.5.1 单元信息查询

API	函数签名	说明
mesh::getElementTypes	(elementTypes)	获取网格中所有 MSH 单元类型标识列表
mesh::getElementProperties	(elementType, elemName, dim, order, numNodes, localNodeCoord, numPrimNodes)	获取单元类型的名称、维度、阶数、节点数和参考单元节点局部坐标
mesh::getElementType	(familyName, order) -> int	根据单元系列名（如 "tetrahedron"）和阶数返回 MSH 类型标识
mesh::getElementsByType	(elementType, elementTags, nodeTags, tag)	按类型获取单元标签和节点（tag<0 获取所有实体）
mesh::getMaxElementTag	(maxTag)	获取当前网格中最大单元标签值

26.5.2 积分点

API	函数签名	说明
mesh::getIntegrationPoints	(elementType, integrationType, localCoord, weights)	获取参考单元中积分点的局部坐标和权重。integrationType 格式为 "家族名+阶数"（如 "Gauss4", "CompositeGauss3"）

26.5.3 雅可比矩阵

API	函数签名	说明
mesh::getJacobians	(elementType, localCoord, jacobians, determinants, coord, tag, task, numTasks)	批量计算所有单元的雅可比矩阵（按列存储）、行列式及积分点物理坐标。支持并行分片
mesh::getJacobian	(elementTag, localCoord, jacobians, determinants, coord)	计算指定单个单元（通过标签）的雅可比数据
mesh::preallocateJacobians	(elementType, numEvalPoints, allocJ, allocDet, allocCoord, jac, det, coord, tag)	预分配雅可比数据缓冲区（多线程并行前必须调用）

26.5.4 基函数

API	函数签名	说明
mesh::getBasisFunctions	(elementType, localCoord, functionSpaceType, numComponents, basisFunctions, numOrientations, wantedOrientations)	获取指定函数空间在积分点处的基函数值或梯度。支持 Lagrange / GradLagrange / H1Legendre / HcurlLegendre 等系列
mesh::getBasisFunctionsOrientation	(elementType, functionSpaceType, basisFunctionsOrientation, tag, task, numTasks)	获取每个单元基函数的方向索引（Lagrange 恒为 0）
mesh::getNumberOfOrientations	(elementType, functionSpaceType) -> int	获取指定单元类型和函数空间的可能方向总数

26.5.5 网格边/面

API	函数签名	说明
mesh::getElementEdgeNodes	(elementType, nodeTags, tag, primary, task, numTasks)	获取单元类型对应所有单元的边节点列表（按单元顺序排列）
mesh::getElementFaceNodes	(elementType, faceType, nodeTags, tag, primary, task, numTasks)	获取单元类型对应所有单元的面节点列表。faceType: 3=三角形, 4=四边形
mesh::createEdges	(dimTags)	为指定实体（空=全部）创建全局唯一边并分配标签
mesh::createFaces	(dimTags)	为指定实体创建全局唯一面并分配标签
mesh::getEdges	(nodeTags, edgeTags, edgeOrientations)	通过节点对列表查询边的全局标签和方向
mesh::getFaces	(faceType, nodeTags, faceTags, faceOrientations)	通过节点列表查询面的全局标签和方向
mesh::getAllEdges	(edgeTags, edgeNodes)	一次性获取网格中所有唯一边的标签和节点
mesh::getAllFaces	(faceType, faceTags, faceNodes)	一次性获取网格中所有唯一面的标签和节点

26.5.6 离散实体与单元添加

API	函数签名	说明
model::addDiscreteEntity	(dim, tag, boundary) -> int	添加一个离散实体（无参数化几何，仅由网格定义），返回实体标签
mesh::addElementsByType	(entityTag, elementType, elementTags, nodeTags)	向指定实体添加一批指定类型的单元

26.5.7 通用选项

选项	说明
mesh::setOrder(N)	设置单元阶数（1=线性, 2=二次,...），影响节点数但不影响几何精度
mesh::setTransfiniteAutomatic()	对所有可结构化剖分的面/体自动应用 transfinite 算法

---

26.6 注意事项与最佳实践

1. 雅可比按列存储：jacobians 向量按列排列 [Jxu, Jyu, Jzu, Jxv, Jyv, Jzv, Jxw, Jyw, Jzw]，而非按行。在构造物理空间梯度 ∂N/∂x = J^{-T} * ∂N/∂u 时务必注意索引的正确转换。

2. 多线程并行：getJacobians() 和 getBasisFunctionsOrientation() 支持 task/numTasks 分片并行。使用前必须调用对应的 preallocate 函数预分配输出缓冲区。

3. 积分规则选择：对三角形/四面体单元，"Gauss" 系列的积分点可能落在参考单元外部（尤其高阶时）。如需保证所有积分点均在单元内部，使用 "CompositeGauss" 张量积规则。

4. 边/面节点顺序：getElementEdgeNodes() 和 getElementFaceNodes() 的返回顺序与 getElements() / getElementsByType() 的单元顺序严格一致。利用这一性质可直接通过整数除法 i / (每单元边数或面数) 建立边/面到单元的索引映射。

5. 全局标签与唯一性：createEdges() / createFaces() 必须在查询边/面标签前调用，否则 getEdges() / getFaces() 无法返回有效标签。getAllEdges() / getAllFaces() 同样依赖此前提。

6. DG 邻接关系：faces2elems 映射中共享计数 = 2 的面是内部面，= 1 的是边界面。在 DG 方法中内部面需计算双侧数值通量，边界面仅需单侧通量（结合边界条件）。

7. 离散实体单元标签：在离散实体上创建面单元时，为避免标签冲突，建议基于 getMaxElementTag() 偏移生成新标签。Gmsh 不允许单元标签重复。

8. 方向与基函数：对于基于边的有限元（如 Nedelec 单元），getEdges() 返回的 edgeOrientations 和 getBasisFunctionsOrientation() 返回的方向索引共同决定基函数在物理单元中的实际朝向，这对 Maxwell 方程的数值求解至关重要。

---

26.7 本章小结

• 本章覆盖了 CAE 前处理中有限元核心数据的获取：积分点（参考单元中的数值积分位置和权重）、雅可比矩阵（参考空间到物理空间的映射和体积变形）、基函数（在积分点处的值和梯度），以及内部边/面（用于 DG 方法、内聚力单元等）。
• getIntegrationPoints() + getJacobians() + getBasisFunctions() 三个 API 组合即可完成标准有限元刚度矩阵组装的全部数据准备。
• getBasisFunctions() 支持 Lagrange、H1Legendre、HcurlLegendre 等多种函数空间，覆盖从标准连续 Galerkin（CG）到 DG 再到 Maxwell 方程求解器的需求。
• createEdges() / createFaces() + getElementEdgeNodes() / getElementFaceNodes() + addElementsByType() 构成了 DG 方法内部面积分的完整前处理流水线。
• 雅可比矩阵的按列存储、边/面的方向标识、多线程预分配等细节在实际工程代码中需要特别关注，避免因为数据布局误解导致刚度矩阵计算错误。

---