CCG 06-01:实验准备、网格构造与拓扑预检

本页只解决一个问题

怎样构造可靠输入,并在求解前拒绝错误拓扑?

本页是原长篇讲义的单概念版本。第一次阅读只追踪本页定义的对象;看见未掌握的前置符号时,先沿链接退回,不要跳过定义硬背结论。

开始前

1. 实验假设

输入网格必须满足:

  • 单连通圆盘拓扑;
  • 恰好一条边界环;
  • 所有面朝向一致;
  • 没有重复面、孤立顶点、非流形边;
  • 每个三角形面积严格为正。

如果输入是球面或环面,先切开;如果有多个洞,需要先设计切图,不能把所有边界随意塞到同一个圆上。

2. 依赖与运行方式

python -m pip install numpy scipy matplotlib

本文代码按出现顺序放在同一个 parameterization.py 中即可运行。

3. 构造一个弯曲圆盘网格

我们从规则矩形三角剖分开始,再将内部顶点弯到三维。原始网格拓扑仍是圆盘,但三维内在边长不再等于平面矩形的边长。

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
from collections import Counter, defaultdict, deque
from scipy.sparse import coo_matrix
from scipy.sparse.linalg import spsolve


def curved_rectangle_mesh(nx=24, ny=16):
    vertices = []
    for j in range(ny + 1):
        y = j / ny
        for i in range(nx + 1):
            x = i / nx
            z = 0.28 * np.sin(np.pi * x) * np.sin(np.pi * y)
            vertices.append((x, y, z))

    def vid(i, j):
        return j * (nx + 1) + i

    faces = []
    for j in range(ny):
        for i in range(nx):
            a, b = vid(i, j), vid(i + 1, j)
            c, d = vid(i, j + 1), vid(i + 1, j + 1)
            # 交错对角线,两个三角形都保持一致方向。
            if (i + j) % 2 == 0:
                faces += [(a, b, d), (a, d, c)]
            else:
                faces += [(a, b, c), (b, d, c)]

    return np.asarray(vertices, float), np.asarray(faces, int)

4. 拓扑预检与边界环

4.1 从有向面恢复边

对一致定向曲面,内部无向边出现两次且方向相反;边界边只出现一次。记录这唯一方向,就得到边界的诱导方向。

def mesh_edges_and_boundary_loops(vertex_count, faces):
    use_count = Counter()
    directed_occurrence = defaultdict(list)

    for face in faces:
        i, j, k = map(int, face)
        for a, b in [(i, j), (j, k), (k, i)]:
            edge = tuple(sorted((a, b)))
            use_count[edge] += 1
            directed_occurrence[edge].append((a, b))

    non_manifold = [e for e, count in use_count.items() if count > 2]
    if non_manifold:
        raise ValueError(f"non-manifold edges: {non_manifold[:8]}")

    inconsistent = []
    for edge, occurrences in directed_occurrence.items():
        if len(occurrences) == 2:
            (a, b), (c, d) = occurrences
            if (a, b) != (d, c):
                inconsistent.append(edge)
    if inconsistent:
        raise ValueError(f"inconsistent face orientation: {inconsistent[:8]}")

    boundary_halfedges = [
        directed_occurrence[e][0]
        for e, count in use_count.items()
        if count == 1
    ]

    outgoing = defaultdict(list)
    incoming = defaultdict(list)
    for a, b in boundary_halfedges:
        outgoing[a].append(b)
        incoming[b].append(a)

    boundary_vertices = set(outgoing) | set(incoming)
    for v in boundary_vertices:
        if len(outgoing[v]) != 1 or len(incoming[v]) != 1:
            raise ValueError(f"boundary is not a 1-manifold at vertex {v}")

    loops = []
    unused = set(boundary_halfedges)
    while unused:
        start_edge = min(unused)
        start, current = start_edge
        loop = [start]
        unused.remove(start_edge)

        while current != start:
            loop.append(current)
            nxt = outgoing[current][0]
            edge = (current, nxt)
            if edge not in unused:
                raise ValueError("boundary loop closes incorrectly")
            unused.remove(edge)
            current = nxt
        loops.append(np.asarray(loop, dtype=int))

    V, E, F = vertex_count, len(use_count), len(faces)
    report = {
        "V": V,
        "E": E,
        "F": F,
        "chi": V - E + F,
        "boundary_components": len(loops),
    }
    return sorted(use_count), loops, report


def validate_disk(vertices, faces):
    edges, loops, report = mesh_edges_and_boundary_loops(len(vertices), faces)
    if report["chi"] != 1 or report["boundary_components"] != 1:
        raise ValueError(f"expected disk topology, got {report}")

    adjacency = defaultdict(list)
    for i, j in edges:
        adjacency[i].append(j)
        adjacency[j].append(i)
    visited = {0}
    queue = deque([0])
    while queue:
        i = queue.popleft()
        for j in adjacency[i]:
            if j not in visited:
                visited.add(j)
                queue.append(j)
    if len(visited) != len(vertices):
        raise ValueError("mesh is disconnected or contains isolated vertices")

    double_areas = []
    for i, j, k in faces:
        area2 = np.linalg.norm(
            np.cross(vertices[j] - vertices[i], vertices[k] - vertices[i])
        )
        double_areas.append(area2)
    if min(double_areas) <= 1e-14:
        raise ValueError("mesh contains a degenerate triangle")
    return edges, loops[0], report

4.2 为什么χ=1,b=1\chi=1,b=1 还不完全够

对连通、可定向流形,

χ=2−2g−b.\chi=2-2g-b.

χ=1,b=1\chi=1,b=1 推出g=0g=0。但如果网格不是流形或不连通,这个结论不能直接使用。因此实际预检还应做连通分量、顶点邻域和面朝向检查;上面代码展示的是本实验需要的最低版本。

读完必须能回答

  1. 不看公式,用自己的话回答开头的核心问题。
  2. 指出本页最重要公式中每个符号的类型:标量、向量、矩阵、映射还是集合。
  3. 找出定理成立所需的假设,并举一个删除某条假设后可能失败的例子。
  4. 对最小例子手算一遍;若本页含代码,再比较手算与程序输出。

前后课程

返回模块目录 | 下一课:边界映射与调和参数化