Skip to content

第12章 同调与上同调计算

12.1 学习目标

  • 理解同调(Homology)和上同调(Cohomology)的直观含义：描述区域中"绕洞的闭合路径"和"切断裂隙的割面"
  • 理解电磁场有限元中"割面(cut)"的作用：将非单连通域切开为单连通域，使标量势可解
  • 掌握 addHomologyRequest() + mesh::generate() 自动计算同调/上同调空间基的方法
  • 学会使用 computeHomology() 在已有网格上触发计算并获取生成的物理组
  • 能解释计算结果：每个生成元对应一个物理组，代表一条割面或闭合环链

12.2 核心概念说明

12.2.1 同调与上同调的直观理解

在有限元前处理中，同调(Homology)和上同调(Cohomology)用于回答：区域内有多少个独立的"洞"或"环"，以及它们如何影响场的求解。

  • 1维同调群 H1：描述空间中"绕洞的独立闭合环路"。一个甜甜圈有两条独立闭合环
  • 2维同调群 H2：描述空间中"包围空腔的独立闭合曲面"
  • 1维上同调群 H1：同调群的对偶，对应"切断所有绕洞环路的割面"(cutting surfaces)
  • 相对同调(Relative Homology)：将某些子区域（如终端面）视为"模掉"对象，在商空间上寻找基

Gmsh 将计算结果存储为物理组(Physical Group)：每个生成元(generator)对应一个物理组，其中的网格单元构成该生成元的代表链(representative chain)。

12.2.2 电磁场问题中的"割面"需求

在电磁场有限元仿真中，经常遇到含孔洞的导体区域。该区域在拓扑上是非单连通(multiply connected)的——存在绕孔的闭合路径使得标量磁势(magnetic scalar potential)产生多值性问题（绕孔一周后势值不闭合）。

解决方法：在区域内插入一个割面(cut)——一个穿过孔洞的假想曲面，将多连通域"切开"为单连通域。在割面两侧施加磁动势跳变条件后，标量势即可唯一确定。

Gmsh 的同调/上同调计算可以自动找出割面的位置并生成对应物理组：

  • 同调基(Homology basis) 给出"薄割面"(thin cut)：沿边界面分布的 2D 单元带
  • 上同调基(Cohomology basis) 给出"厚割面"(thick cut)：较宽的 2D 带状区域，更适合通用有限元求解

12.2.3 示例几何概述

本章示例：一个 10x10x2 的长方体导体，中心含一个 2x2 的方形贯穿孔洞，左右两侧各有 2 个终端面（模拟电流入口/出口）。该区域有 1 个独立绕孔回路，需要 1 个割面。

12.3 C++ 代码逐段讲解

12.3.1 初始化与几何点

#include <set>
#include <algorithm>
#include <gmsh.h>

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

  double m = 0.5;  // 网格尺寸
  double h = 2.0;  // z 方向拉伸高度

  // 外边界四角点 (10x10)
  gmsh::model::geo::addPoint(0, 0, 0, m, 1);
  gmsh::model::geo::addPoint(10, 0, 0, m, 2);
  gmsh::model::geo::addPoint(10, 10, 0, m, 3);
  gmsh::model::geo::addPoint(0, 10, 0, m, 4);

  // 内部方形孔洞角点 (2x2, 中心5,5)
  gmsh::model::geo::addPoint(4, 4, 0, m, 5);
  gmsh::model::geo::addPoint(6, 4, 0, m, 6);
  gmsh::model::geo::addPoint(6, 6, 0, m, 7);
  gmsh::model::geo::addPoint(4, 6, 0, m, 8);

  // 终端分界点(将左右边各分为两个终端面)
  gmsh::model::geo::addPoint(2, 0, 0, m, 9);
  gmsh::model::geo::addPoint(8, 0, 0, m, 10);
  gmsh::model::geo::addPoint(2, 10, 0, m, 11);
  gmsh::model::geo::addPoint(8, 10, 0, m, 12);

addPoint(x, y, z, meshSize, tag) 创建几何点。前三个参数为坐标，第四个参数为目标网格尺寸，第五个参数为点标签（可自定义）。点 5-8 定义内部孔洞，点 9-12 位于外边界上用于切分终端面。

12.3.2 线段、曲线环与平面

  // 外边界线段 (按逆时针方向)
  gmsh::model::geo::addLine(1, 9, 1);    gmsh::model::geo::addLine(9, 10, 2);
  gmsh::model::geo::addLine(10, 2, 3);   gmsh::model::geo::addLine(2, 3, 4);
  gmsh::model::geo::addLine(3, 12, 5);   gmsh::model::geo::addLine(12, 11, 6);
  gmsh::model::geo::addLine(11, 4, 7);   gmsh::model::geo::addLine(4, 1, 8);

  // 内部孔洞线段 (顺时针方向)
  gmsh::model::geo::addLine(5, 6, 9);    gmsh::model::geo::addLine(6, 7, 10);
  gmsh::model::geo::addLine(7, 8, 11);   gmsh::model::geo::addLine(8, 5, 12);

  // 曲线环: 外环13(逆时针) + 内环14(顺时针)
  gmsh::model::geo::addCurveLoop({6, 7, 8, 1, 2, 3, 4, 5}, 13);
  gmsh::model::geo::addCurveLoop({11, 12, 9, 10}, 14);
  // 带孔洞的平面: 外环在前, 内环在后
  gmsh::model::geo::addPlaneSurface({13, 14}, 15);

addCurveLoop 将线段组成闭合环。addPlaneSurface 中第一个环为外边界（逆时针），第二个环为孔洞（顺时针），Gmsh 会自动处理方向。这样就构建了一个带 2x2 方孔的 10x10 平面。

12.3.3 拉伸生成三维体

  std::vector<std::pair<int, int>> e;
  gmsh::model::geo::extrude({{2, 15}}, 0, 0, h, e);
  gmsh::model::geo::synchronize();

extrude 将平面沿 z 轴拉伸 h=2，生成三维体及其 6 个侧面。返回值 e 中的关键条目：

索引 内容 用途
e[1] 3D 体 计算域
e[3] 左侧下终端面 终端物理组
e[5] 左侧上终端面 终端物理组
e[7] 右侧上终端面 终端物理组
e[9] 右侧下终端面 终端物理组

synchronize() 将几何内核数据同步到 model 模块，后续查询操作才能正确执行。

12.3.4 定义物理组

物理组是同调计算的基石，用于指定计算域和相对子域。

  // 计算域: 三维体
  int domain_physical_tag = 1001;
  gmsh::model::addPhysicalGroup(3, {e[1].second}, domain_physical_tag,
                                "Whole domain");

  // 四个终端面: 相对同调的子域
  std::vector<int> terminal_tags = {e[3].second, e[5].second,
                                    e[7].second, e[9].second};
  int terminals_physical_tag = 2001;
  gmsh::model::addPhysicalGroup(2, terminal_tags, terminals_physical_tag,
                                "Terminals");

  // 获取域的全部边界面，构造"补集"(Complement = 边界面 - 终端面)
  std::vector<std::pair<int, int>> boundary_dimtags;
  gmsh::model::getBoundary({{3, e[1].second}}, boundary_dimtags, false, false);

  std::vector<int> boundary_tags, complement_tags;
  for (auto b : boundary_dimtags) {
    boundary_tags.push_back(b.second);
    complement_tags.push_back(b.second);
  }
  for (auto t : terminal_tags) {
    auto it = std::find(complement_tags.begin(), complement_tags.end(), t);
    if (it != complement_tags.end()) complement_tags.erase(it);
  }

  int complement_physical_tag = 2003;
  gmsh::model::addPhysicalGroup(2, complement_tags, complement_physical_tag,
                                "Complement");

addPhysicalGroup(dim, tags, tag, name) 将实体归入命名物理组。getBoundary 获取 3D 体的包围面。建立三个关键物理组：

物理组 维度 说明
Whole domain (1001) 3D 同调计算域
Terminals (2001) 2D 相对同调的子域，模掉对象
Complement (2003) 2D 除终端面外的所有外表面

12.3.5 提交同调/上同调计算请求

  // 相对同调: 域模掉终端面 → 薄割面(thin cut)
  gmsh::model::mesh::addHomologyRequest(
      "Homology",
      {domain_physical_tag},
      {terminals_physical_tag},
      {0, 1, 2, 3});

  // 相对同调: 域模掉补集 → 同构的薄割面(位于非终端外表面)
  gmsh::model::mesh::addHomologyRequest(
      "Homology",
      {domain_physical_tag},
      {complement_physical_tag},
      {0, 1, 2, 3});

  // 相对上同调: 域模掉终端面 → 厚割面(thick cut)
  gmsh::model::mesh::addHomologyRequest(
      "Cohomology",
      {domain_physical_tag},
      {terminals_physical_tag},
      {0, 1, 2, 3});

addHomologyRequest(type, domainTags, subdomainTags, dims) 参数详解：

参数 含义 默认值
type "Homology"=同调基, "Cohomology"=上同调基 "Homology"
domainTags 计算域物理组标签列表；空=整个网格 {}
subdomainTags 相对同调子域标签列表；空=绝对同调 {}
dims 感兴趣的维度列表；空=全部维度(0-3) {}

三次调用的意义： - 第 1 次：计算相对同调 H(域, 终端)，基链位于终端面上 - 第 2 次：计算 H(域, 补集)，基链位于非终端外表面上（同构但表示不同） - 第 3 次：计算相对上同调 H(域, 终端)，生成"厚割面"，适合有限元求解器使用

12.3.6 网格生成与计算执行

  gmsh::model::mesh::generate(3);

mesh::generate(3) 生成三维网格，并在流水线末端自动执行所有 addHomologyRequest() 提交的计算。计算完成后，Gmsh 为同调/上同调空间的每个生成元自动创建新的物理组并填充代表链单元。

如果需要在网格生成后单独触发计算（例如对已有网格），可以使用：

  std::vector<std::pair<int, int>> newPhysicals;
  gmsh::model::mesh::computeHomology(newPhysicals);
  // newPhysicals 即自动创建的物理组 (dim, tag) 列表

computeHomology() 执行所有待处理的(共)同调请求，并返回新生成的物理组标签。注意 Gmsh C++ API 中并没有单独的 computeCohomology() 函数——上同调通过 addHomologyRequest("Cohomology", ...) 指定类型，由 computeHomology() 统一执行。

12.3.7 输出与可视化

  gmsh::write("t14.msh");

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

  gmsh::finalize();
  return 0;
}

write("t14.msh") 将网格及所有物理组（含自动生成的同调/上同调物理组）写入文件。GUI 中可在左侧 Visibility 面板查看自动生成的物理组——高亮后即可看到割面的空间位置与形状。

结果的物理组命名格式： - 绝对同调："Homology: DomainTag_dim" - 相对同调："Homology: DomainTag|SubdomainTag_dim" - 相对上同调："Cohomology: DomainTag|SubdomainTag_dim"

其中 dim 为 0/1/2/3，对应各维度的生成元。

12.4 完整可运行代码

// =============================================================================
// 第12章 同调与上同调计算 — 完整可运行代码
//
// 编译: g++ -std=c++11 -o ch12 ch12.cpp -lgmsh
// 运行: ./ch12          (GUI 查看结果)
//       ./ch12 -nopopup (仅输出 t14.msh)
// =============================================================================

#include <set>
#include <algorithm>
#include <gmsh.h>

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

  double m = 0.5, h = 2.0;

  // -- 点 --
  gmsh::model::geo::addPoint(0, 0, 0, m, 1);
  gmsh::model::geo::addPoint(10, 0, 0, m, 2);
  gmsh::model::geo::addPoint(10, 10, 0, m, 3);
  gmsh::model::geo::addPoint(0, 10, 0, m, 4);
  gmsh::model::geo::addPoint(4, 4, 0, m, 5);
  gmsh::model::geo::addPoint(6, 4, 0, m, 6);
  gmsh::model::geo::addPoint(6, 6, 0, m, 7);
  gmsh::model::geo::addPoint(4, 6, 0, m, 8);
  gmsh::model::geo::addPoint(2, 0, 0, m, 9);
  gmsh::model::geo::addPoint(8, 0, 0, m, 10);
  gmsh::model::geo::addPoint(2, 10, 0, m, 11);
  gmsh::model::geo::addPoint(8, 10, 0, m, 12);

  // -- 线 --
  gmsh::model::geo::addLine(1, 9, 1);   gmsh::model::geo::addLine(9, 10, 2);
  gmsh::model::geo::addLine(10, 2, 3);  gmsh::model::geo::addLine(2, 3, 4);
  gmsh::model::geo::addLine(3, 12, 5);  gmsh::model::geo::addLine(12, 11, 6);
  gmsh::model::geo::addLine(11, 4, 7);  gmsh::model::geo::addLine(4, 1, 8);
  gmsh::model::geo::addLine(5, 6, 9);   gmsh::model::geo::addLine(6, 7, 10);
  gmsh::model::geo::addLine(7, 8, 11);  gmsh::model::geo::addLine(8, 5, 12);

  // -- 曲线环与平面 --
  gmsh::model::geo::addCurveLoop({6, 7, 8, 1, 2, 3, 4, 5}, 13);
  gmsh::model::geo::addCurveLoop({11, 12, 9, 10}, 14);
  gmsh::model::geo::addPlaneSurface({13, 14}, 15);

  // -- 拉伸 --
  std::vector<std::pair<int, int>> e;
  gmsh::model::geo::extrude({{2, 15}}, 0, 0, h, e);
  gmsh::model::geo::synchronize();

  // -- 物理组 --
  int domTag = 1001, termTag = 2001, compTag = 2003;
  gmsh::model::addPhysicalGroup(3, {e[1].second}, domTag, "Whole domain");

  std::vector<int> term = {e[3].second, e[5].second, e[7].second, e[9].second};
  gmsh::model::addPhysicalGroup(2, term, termTag, "Terminals");

  std::vector<std::pair<int, int>> bnd;
  gmsh::model::getBoundary({{3, e[1].second}}, bnd, false, false);

  std::vector<int> comp;
  for (auto b : bnd) comp.push_back(b.second);
  for (auto t : term) {
    auto it = std::find(comp.begin(), comp.end(), t);
    if (it != comp.end()) comp.erase(it);
  }
  gmsh::model::addPhysicalGroup(2, comp, compTag, "Complement");

  // -- 同调/上同调请求 --
  // 薄割面: 相对同调 H(域, 终端)
  gmsh::model::mesh::addHomologyRequest("Homology", {domTag}, {termTag},
                                        {0, 1, 2, 3});
  // 薄割面: 相对同调 H(域, 补集)
  gmsh::model::mesh::addHomologyRequest("Homology", {domTag}, {compTag},
                                        {0, 1, 2, 3});
  // 厚割面: 相对上同调 H(域, 终端)
  gmsh::model::mesh::addHomologyRequest("Cohomology", {domTag}, {termTag},
                                        {0, 1, 2, 3});

  // -- 网格生成（自动执行同调计算） --
  gmsh::model::mesh::generate(3);

  gmsh::write("t14.msh");

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

  gmsh::finalize();
  return 0;
}

12.5 关键 API 速查表

API 说明
gmsh::model::mesh::addHomologyRequest(type, domainTags, subdomainTags, dims) 提交同调/上同调计算请求。若在 mesh::generate() 之前调用，计算在网格划分尾部自动执行。type="Homology" 或 "Cohomology"
gmsh::model::mesh::computeHomology(dimTags) 执行所有通过 addHomologyRequest() 提交的(共)同调计算。返回 vectorpair 为新创建的物理组 (dim, tag) 列表。可在网格生成后单独调用
gmsh::model::mesh::clearHomologyRequests() 清除所有待处理的同调/上同调计算请求
gmsh::model::addPhysicalGroup(dim, tags, tag, name) 创建物理组。dim 为维度(1 边, 2 面, 3 体)，tags 为实体标签列表，tag 为物理组标签，name 为可选名称
gmsh::model::getBoundary(dimTags, outDimTags, combined, oriented) 获取实体边界。combined=true 合并连续边界，oriented=true 保留方向信息
gmsh::model::geo::addPoint(x, y, z, meshSize, tag) 创建几何点。前三个参数为坐标，meshSize 为目标网格尺寸
gmsh::model::geo::addLine(startTag, endTag, tag) 在两点间创建直线段
gmsh::model::geo::addCurveLoop(lineTags, tag) 由线段列表创建闭合曲线环
gmsh::model::geo::addPlaneSurface(loopTags, tag) 由曲线环创建平面。第一个环为外边界（逆时针），后续环为孔洞（顺时针）
gmsh::model::geo::extrude(dimTags, dx, dy, dz, outDimTags) 沿 (dx, dy, dz) 拉伸实体，返回新生成实体列表
gmsh::model::geo::synchronize() 将几何内核同步到 model 模块（网格生成/查询前的必要步骤）
gmsh::model::mesh::generate(dim) 生成 dim 维网格，并自动执行已提交的同调计算请求

参考文献：M. Pellikka, S. Suuriniemi, L. Kettunen and C. Geuzaine. Homology and Cohomology Computation in Finite Element Modeling. SIAM Journal on Scientific Computing, 35(5), pp. 1195-1214, 2013.