9.3 四面体重网格化 (Tetrahedral_remeshing)

上一节: 9.2 曲面网格生成

9.3.1 数学理论基础

重网格化问题定义

四面体重网格化（Tetrahedral Remeshing）是改善现有四面体网格质量的过程。给定一个四面体网格 $M$ ，重网格化的目标是：

改善质量：提高四面体的形状质量
调整尺寸：根据尺寸场调整网格密度
保持特征：保持几何特征（边、角）
保持拓扑：保持原始网格的拓扑结构

质量度量

1. 半径-边比（Radius-Edge Ratio）

$ρ = \frac{R}{L _{min}}$

其中 $R$ 是外接球半径， $L_{min}$ 是最短边长。

等边四面体： $ρ = \frac{6}{4} \approx 0.612$
退化四面体： $ρ \to \infty$

2. 二面角（Dihedral Angle）

四面体两个面之间的夹角，范围 $(0, π)$ 。

等边四面体： $arccos (1/3) \approx 70.5 3^{\circ}$
薄片（Sliver）：二面角接近 $0^{\circ}$ 或 $18 0^{\circ}$

3. 体积-边比（Volume-Edge Ratio）

$η = \frac{V}{L _{r m s}^{3}}$

其中 $V$ 是体积， $L_{r m s}$ 是边长的均方根。

重网格化操作

CGAL的Tetrahedral_remeshing基于局部网格修改操作：

1. 边分裂（Edge Split）

将一条边在中点处分裂，连接相邻顶点形成新四面体。

    v3                    v3
   /|\                   /|\
  / | \                 / | \
 /  |  \               /  |  \
v1--*--v2   -->    v1--c--v2
 \  |  /               \  |  /
  \ | /                 \ | /
   \|/                   \|/
    v4                    v4

2. 边折叠（Edge Collapse）

将一条边折叠为一个顶点，删除相邻四面体。

3. 边翻转（Edge Flip）

改变边的连接方式，改变局部拓扑。

2-3翻转：两个四面体 -> 三个四面体
3-2翻转：三个四面体 -> 两个四面体（逆操作）
4-4翻转：四个四面体重新排列

4. 顶点重定位（Vertex Relocation）

移动顶点位置以改善局部质量。

Laplacian平滑： $v_{n e w} = \frac{1}{∣ N ( v ) ∣} \sum_{u \in N (v)} u$

优化平滑：求解局部优化问题

尺寸场自适应

重网格化支持自适应尺寸场 $ψ (x)$ ：

分裂过大边：当边长 $> ψ$ 时分裂
折叠过小边：当边长 $< ψ / 2$ 时折叠
保持目标尺寸：边长维持在 $[ψ / 2, ψ]$ 范围内

9.3.2 架构分析

整体架构

Tetrahedral_remeshing
├── Remeshing_triangulation（重网格化三角剖分）
│   ├── Tetrahedral_remeshing_triangulation
│   └── Remeshing_vertex_base
│
├── Remeshing_algorithm（重网格化算法）
│   ├── Split_edges（边分裂）
│   ├── Collapse_edges（边折叠）
│   ├── Flip_edges（边翻转）
│   └── Move_vertices（顶点移动）
│
├── Sizing_field（尺寸场）
│   ├── Constant_sizing_field
│   ├── Adaptive_sizing_field
│   └── Custom_sizing_field
│
└── Criteria（质量标准）
    ├── Dihedral_angle_bound
    └── Cell_quality_threshold

核心类

1. 重网格化三角剖分

template <class Gt, class Vb = Default, class Cb = Default>
class Tetrahedral_remeshing_triangulation
    : public Delaunay_triangulation_3<Gt, Tds> {
public:
    // 额外的重网格化数据
    bool is_in_domain(Cell_handle c) const;
    void set_subdomain(Cell_handle c, Subdomain_index idx);
    
    // 边访问
    Edge_iterator edges_in_domain_begin() const;
    Edge_iterator edges_in_domain_end() const;
};

2. 重网格化算法

template <class Tr>
class Tetrahedral_remeshing {
public:
    // 构造函数
    Tetrahedral_remeshing(Tr& tr, 
                          const Sizing_field& sizing,
                          const Remeshing_options& options);
    
    // 执行重网格化
    void remesh();
    
    // 分步执行
    void split_edges();
    void collapse_edges();
    void flip_edges();
    void smooth_vertices();
};

3. 尺寸场概念

class Sizing_field {
public:
    // 必需的操作
    FT operator()(const Point_3& p) const;  // 在点p处的目标尺寸
    FT operator()(const Edge& e) const;      // 在边e上的目标尺寸
};

算法流程

输入：四面体网格 M，尺寸场 ψ，迭代次数 N
输出：改善后的网格 M'

for iter = 1 to N:
    // 1. 边分裂
    for each edge e in M:
        if length(e) > ψ(e):
            split_edge(e)
    
    // 2. 边折叠
    for each edge e in M:
        if length(e) < ψ(e)/√2:
            if collapse_preserves_topology(e):
                collapse_edge(e)
    
    // 3. 边翻转
    for each edge e in M:
        if flip_improves_quality(e):
            flip_edge(e)
    
    // 4. 顶点平滑
    for each vertex v in M:
        if v is not on feature:
            v_new = optimize_position(v)
            if quality_improves(v_new):
                move_vertex(v, v_new)

// 5. 投影到边界
project_boundary_vertices_to_surface()

9.3.3 实现细节

1. 初始化三角剖分

#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation.h>
 
// 类型定义
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation<K> Remeshing_triangulation;
 
typedef Remeshing_triangulation::Point Point;
typedef Remeshing_triangulation::Vertex_handle Vertex_handle;
typedef Remeshing_triangulation::Cell_handle Cell_handle;
 
// 创建三角剖分
Remeshing_triangulation tr;
 
// 从现有网格初始化
// 方法1：从顶点-单元列表
std::vector<Point> points;
std::vector<std::array<int, 4>> cells;
// ... 填充点和单元 ...
 
for (const auto& p : points) {
    tr.insert(p);
}
 
// 标记域
for (auto cit = tr.finite_cells_begin(); cit != tr.finite_cells_end(); ++cit) {
    tr.set_subdomain(cit, 1);  // 标记为域内
}

2. 定义尺寸场

#include <CGAL/Tetrahedral_remeshing/Sizing_field.h>
 
// 常数尺寸场
struct Constant_sizing {
    double target_size;
    
    double operator()(const Point& p) const {
        return target_size;
    }
    
    double operator()(const Remeshing_triangulation::Edge& e) const {
        return target_size;
    }
};
 
// 自适应尺寸场
struct Adaptive_sizing {
    double operator()(const Point& p) const {
        double r = std::sqrt(p.x()*p.x() + p.y()*p.y() + p.z()*p.z());
        return 0.1 / (1.0 + r);  // 远离原点处更粗
    }
};
 
// 从背景网格插值
typedef CGAL::Tetrahedral_remeshing::Interpolated_sizing_field<K> Sizing_field;
Sizing_field sizing(background_mesh);

3. 执行重网格化

#include <CGAL/tetrahedral_remeshing.h>
 
// 设置选项
CGAL::Tetrahedral_remeshing::Remeshing_options options;
options.number_of_iterations = 3;
options.preserve_features = true;
options.smooth_constrained_edges = false;
 
// 执行重网格化
Constant_sizing sizing(0.1);  // 目标尺寸0.1
CGAL::tetrahedral_remeshing(tr, sizing, options);

4. 特征保护

// 标记特征边
void mark_features(Remeshing_triangulation& tr) {
    // 基于二面角检测特征边
    for (auto eit = tr.finite_edges_begin(); eit != tr.finite_edges_end(); ++eit) {
        double dihedral_angle = compute_dihedral_angle(tr, *eit);
        
        if (dihedral_angle < CGAL::to_radians(120.0)) {
            // 标记为特征边
            tr.set_feature(*eit, true);
        }
    }
}
 
// 标记边界边
void mark_boundary(Remeshing_triangulation& tr) {
    for (auto eit = tr.finite_edges_begin(); eit != tr.finite_edges_end(); ++eit) {
        if (tr.is_boundary(*eit)) {
            tr.set_feature(*eit, true);
        }
    }
}

5. 质量评估

#include <CGAL/Tetrahedral_remeshing/Quality.h>
 
// 计算四面体质量
template <class Cell_handle>
double compute_radius_edge_ratio(Cell_handle c) {
    // 外接球半径
    auto circumsphere = circumcenter(c);
    double R = std::sqrt(CGAL::squared_distance(
        circumsphere.first, c->vertex(0)->point()));
    
    // 最短边
    double L_min = std::numeric_limits<double>::max();
    for (int i = 0; i < 4; i++) {
        for (int j = i+1; j < 4; j++) {
            double L = std::sqrt(CGAL::squared_distance(
                c->vertex(i)->point(), c->vertex(j)->point()));
            L_min = std::min(L_min, L);
        }
    }
    
    return R / L_min;
}
 
// 统计网格质量
void analyze_quality(const Remeshing_triangulation& tr) {
    double min_ratio = std::numeric_limits<double>::max();
    double max_ratio = 0;
    double sum_ratio = 0;
    int count = 0;
    
    for (auto cit = tr.finite_cells_begin(); cit != tr.finite_cells_end(); ++cit) {
        if (tr.is_in_domain(cit)) {
            double ratio = compute_radius_edge_ratio(cit);
            min_ratio = std::min(min_ratio, ratio);
            max_ratio = std::max(max_ratio, ratio);
            sum_ratio += ratio;
            count++;
        }
    }
    
    std::cout << "质量统计:" << std::endl;
    std::cout << "  最小半径-边比: " << min_ratio << std::endl;
    std::cout << "  最大半径-边比: " << max_ratio << std::endl;
    std::cout << "  平均半径-边比: " << sum_ratio / count << std::endl;
}

9.3.4 代码示例

完整示例：网格优化

#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Tetrahedral_remeshing/Remeshing_triangulation.h>
#include <CGAL/tetrahedral_remeshing.h>
#include <CGAL/IO/File_medit.h>
#include <iostream>
#include <fstream>
 
// 类型定义
typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
typedef CGAL::Tetrahedral_remeshing::Remeshing_triangulation<K> Tr;
 
typedef Tr::Point Point;
typedef Tr::Vertex_handle Vertex_handle;
typedef Tr::Cell_handle Cell_handle;
 
// 尺寸场
struct Sizing_field {
    double target_size;
    
    Sizing_field(double size) : target_size(size) {}
    
    double operator()(const Point& p) const {
        return target_size;
    }
    
    double operator()(const Tr::Edge& e) const {
        return target_size;
    }
};
 
int main(int argc, char* argv[]) {
    // 读取输入网格
    Tr tr;
    std::ifstream input(argv[1]);
    CGAL::IO::read_MEDIT(input, tr);
    input.close();
    
    std::cout << "输入网格:" << std::endl;
    std::cout << "  顶点数: " << tr.number_of_vertices() << std::endl;
    std::cout << "  单元数: " << tr.number_of_finite_cells() << std::endl;
    
    // 分析初始质量
    analyze_quality(tr);
    
    // 设置尺寸场
    Sizing_field sizing(0.1);  // 目标尺寸0.1
    
    // 设置选项
    CGAL::Tetrahedral_remeshing::Remeshing_options options;
    options.number_of_iterations = 3;
    options.preserve_features = true;
    
    // 执行重网格化
    CGAL::tetrahedral_remeshing(tr, sizing, options);
    
    std::cout << "\n输出网格:" << std::endl;
    std::cout << "  顶点数: " << tr.number_of_vertices() << std::endl;
    std::cout << "  单元数: " << tr.number_of_finite_cells() << std::endl;
    
    // 分析最终质量
    analyze_quality(tr);
    
    // 保存结果
    std::ofstream output("remeshed.mesh");
    CGAL::IO::write_MEDIT(output, tr);
    output.close();
    
    return 0;
}

自适应重网格化

// 基于曲率的自适应尺寸场
struct Curvature_based_sizing {
    Tr& tr;
    double min_size, max_size;
    
    double operator()(const Point& p) const {
        // 估计局部曲率（简化示例）
        double curvature = estimate_curvature(tr, p);
        
        // 高曲率区域使用更细的网格
        double size = max_size / (1.0 + curvature);
        return std::clamp(size, min_size, max_size);
    }
    
    double estimate_curvature(const Tr& tr, const Point& p) const {
        // 查找最近顶点
        Vertex_handle v = tr.nearest_vertex(p);
        
        // 基于邻域法向变化估计曲率
        std::vector<Vector> normals;
        auto vc = tr.incident_cells(v);
        auto done = vc;
        
        do {
            if (!tr.is_infinite(vc)) {
                normals.push_back(compute_cell_normal(vc));
            }
            ++vc;
        } while (vc != done);
        
        // 计算法向变化
        double curvature = 0;
        for (size_t i = 0; i < normals.size(); i++) {
            for (size_t j = i+1; j < normals.size(); j++) {
                curvature = std::max(curvature, 
                    1.0 - std::abs(normals[i] * normals[j]));
            }
        }
        
        return curvature;
    }
};

多域重网格化

// 多域尺寸场
struct Multi_domain_sizing {
    std::map<int, double> subdomain_sizes;
    
    double operator()(const Tr::Cell_handle c) const {
        int subdomain = c->subdomain_index();
        auto it = subdomain_sizes.find(subdomain);
        if (it != subdomain_sizes.end()) {
            return it->second;
        }
        return 0.1;  // 默认尺寸
    }
};
 
void multi_domain_remeshing(Tr& tr) {
    Multi_domain_sizing sizing;
    sizing.subdomain_sizes[1] = 0.05;  // 域1使用细网格
    sizing.subdomain_sizes[2] = 0.2;   // 域2使用粗网格
    
    // 标记域间界面为特征
    for (auto fit = tr.finite_facets_begin(); fit != tr.finite_facets_end(); ++fit) {
        Cell_handle c1 = fit->first;
        Cell_handle c2 = c1->neighbor(fit->second);
        
        if (!tr.is_infinite(c2) && c1->subdomain_index() != c2->subdomain_index()) {
            // 标记界面边为特征
            for (int i = 0; i < 3; i++) {
                int vi = (fit->second + i + 1) % 4;
                int vj = (fit->second + (i+1)%3 + 1) % 4;
                Tr::Edge e(c1, vi, vj);
                tr.set_feature(e, true);
            }
        }
    }
    
    CGAL::tetrahedral_remeshing(tr, sizing);
}

迭代优化

void iterative_optimization(Tr& tr, int max_iterations) {
    double target_quality = 2.0;  // 目标半径-边比
    
    for (int iter = 0; iter < max_iterations; iter++) {
        std::cout << "迭代 " << iter + 1 << std::endl;
        
        // 计算当前质量
        double worst_quality = compute_worst_quality(tr);
        std::cout << "  最差质量: " << worst_quality << std::endl;
        
        if (worst_quality < target_quality) {
            std::cout << "  达到目标质量，停止" << std::endl;
            break;
        }
        
        // 自适应尺寸场
        auto sizing = create_adaptive_sizing(tr, target_quality);
        
        // 执行重网格化
        CGAL::Tetrahedral_remeshing::Remeshing_options options;
        options.number_of_iterations = 1;
        CGAL::tetrahedral_remeshing(tr, sizing, options);
    }
}

9.3.5 复杂度分析

时间复杂度

操作	复杂度	说明
边分裂	$O (1)$	局部操作
边折叠	$O (d)$	d为边度数
边翻转	$O (1)$	局部操作
顶点平滑	$O (d)$	d为顶点度数
单次迭代	$O (n)$	n为网格规模
完整重网格化	$O (k \cdot n)$	k为迭代次数

空间复杂度

三角剖分存储： $O (n)$
特征标记： $O (n)$
临时数据结构： $O (n)$

总空间复杂度： $O (n)$

收敛性

经验观察：

通常3-5次迭代即可达到稳定状态
每次迭代显著改善最差质量单元
平均质量单调改善

9.3.6 理论保证

拓扑保持

重网格化算法保证：

同胚保持：输出网格与输入网格同胚
边界保持：边界拓扑不变
特征保持：标记的特征边/角保持

质量改善

经验保证：

最差半径-边比通常改善50%以上
平均质量显著改善
消除大部分薄片单元

限制：

无法保证达到任意目标质量
受限于输入网格拓扑
特征约束可能限制优化空间

9.3.7 与其他组件的关系

Tetrahedral_remeshing
├── 依赖
│   ├── Delaunay_triangulation_3（基础三角剖分）
│   ├── Mesh_3（网格生成）
│   └── Triangulation_3（三角剖分操作）
│
├── 用于
│   ├── Mesh_3输出优化
│   ├── 有限元预处理
│   └── 多尺度模拟
│
└── 相关工具
    ├── CGAL::Polygon_mesh_processing（表面重网格化）
    └── CGAL::Mesh_3（初始网格生成）

9.3.8 应用

1. 有限元网格优化

// 优化用于FEM的网格
void optimize_fem_mesh(Tr& tr) {
    // 基于应力场的自适应尺寸
    Stress_based_sizing sizing(stress_field);
    
    // 保持几何特征
    mark_geometric_features(tr);
    
    // 重网格化
    CGAL::tetrahedral_remeshing(tr, sizing);
}

2. 多尺度网格

// 创建多尺度网格
void create_multiscale_mesh(Tr& tr) {
    // 近场细网格，远场粗网格
    Distance_from_center_sizing sizing;
    sizing.center = Point(0, 0, 0);
    sizing.near_size = 0.01;
    sizing.far_size = 0.5;
    sizing.transition_distance = 1.0;
    
    CGAL::tetrahedral_remeshing(tr, sizing);
}

3. 网格粗化

// 粗化过细的网格
void coarsen_mesh(Tr& tr, double target_size) {
    Constant_sizing sizing(target_size);
    
    CGAL::Tetrahedral_remeshing::Remeshing_options options;
    options.number_of_iterations = 2;
    
    CGAL::tetrahedral_remeshing(tr, sizing, options);
}

4. 网格细化

// 在特定区域细化网格
void refine_region(Tr& tr, const std::function<bool(Point)>& region) {
    Region_sizing sizing;
    sizing.region = region;
    sizing.in_region_size = 0.01;
    sizing.out_region_size = 0.1;
    
    CGAL::tetrahedral_remeshing(tr, sizing);
}

Chunibyo

Explorer

9.3 四面体重网格化 (Tetrahedral_remeshing)

9.3 四面体重网格化 (Tetrahedral_remeshing) §

9.3.1 数学理论基础 §

重网格化问题定义 §

质量度量 §

1. 半径-边比（Radius-Edge Ratio） §

2. 二面角（Dihedral Angle） §

3. 体积-边比（Volume-Edge Ratio） §

重网格化操作 §

1. 边分裂（Edge Split） §

2. 边折叠（Edge Collapse） §

3. 边翻转（Edge Flip） §

4. 顶点重定位（Vertex Relocation） §

尺寸场自适应 §

9.3.2 架构分析 §

整体架构 §

核心类 §

1. 重网格化三角剖分 §

2. 重网格化算法 §

3. 尺寸场概念 §

算法流程 §

9.3.3 实现细节 §

1. 初始化三角剖分 §

2. 定义尺寸场 §

3. 执行重网格化 §

4. 特征保护 §

5. 质量评估 §

9.3.4 代码示例 §

完整示例：网格优化 §

自适应重网格化 §

多域重网格化 §

迭代优化 §

9.3.5 复杂度分析 §

时间复杂度 §

空间复杂度 §

收敛性 §

9.3.6 理论保证 §

拓扑保持 §

质量改善 §

9.3.7 与其他组件的关系 §

9.3.8 应用 §

1. 有限元网格优化 §

2. 多尺度网格 §

3. 网格粗化 §

4. 网格细化 §

Graph View

Table of Contents

Backlinks

9.3 四面体重网格化 (Tetrahedral_remeshing)

9.3.1 数学理论基础

重网格化问题定义

质量度量

1. 半径-边比（Radius-Edge Ratio）

2. 二面角（Dihedral Angle）

3. 体积-边比（Volume-Edge Ratio）

重网格化操作

1. 边分裂（Edge Split）

2. 边折叠（Edge Collapse）

3. 边翻转（Edge Flip）

4. 顶点重定位（Vertex Relocation）

尺寸场自适应

9.3.2 架构分析

整体架构

核心类

1. 重网格化三角剖分

2. 重网格化算法

3. 尺寸场概念

算法流程

9.3.3 实现细节

1. 初始化三角剖分

2. 定义尺寸场

3. 执行重网格化

4. 特征保护

5. 质量评估

9.3.4 代码示例

完整示例：网格优化

自适应重网格化

多域重网格化

迭代优化

9.3.5 复杂度分析

时间复杂度

空间复杂度

收敛性

9.3.6 理论保证

拓扑保持

质量改善

9.3.7 与其他组件的关系

9.3.8 应用

1. 有限元网格优化

2. 多尺度网格

3. 网格粗化

4. 网格细化