几何创建

按照src目录下的shared、for_2d_build和for_3d_build分为二维与三维通用方法、二维专用方法和三维专用方法分别介绍。最后，还会介绍空shape的使用。

二维与三维通用方法

SPHinXsys自己开发的，不依赖第三方库。

用GeometricShape创建

只能创建Box（二维是矩形，三维是长方体）和Ball（二维是圆，三维是球）。

创建Box有三种方式，参考以下三个构造函数：

// 指定Box中心点、半尺寸和名字。因为默认以原点为中心，所以这里的transform就是平移到指定中心的矢量。
GeometricShapeBox::GeometricShapeBox(const Transform &transform, const Vecd &halfsize,
                      const std::string &name = "GeometricShapeBox");
// 指定整个bounding_box为形状，指定名字
explicit GeometricShapeBox::GeometricShapeBox(const BoundingBox &bounding_box,
                      const std::string &name = "GeometricShapeBox");
// 指定最小坐标点、最大坐标点和名字
GeometricShapeBox::GeometricShapeBox(const Vecd &lower_bound, const Vecd &upper_bound,
                      const std::string &name = "GeometricShapeBox");

创建Ball只有一种方式：

// 指定Ball中心点、半径和名字
explicit GeometricShapeBall::GeometricShapeBall(const Vecd &center, Real radius,
                      const std::string &name = "GeometricShapeBall");

平移与旋转变换（Transform / Rotation）

在很多创建函数里会看到 Transform / translation / “轴线方向”等参数，它们本质上都是在描述几何体如何从“局部坐标系（shape frame）”映射到“全局坐标系（base/world）”。

Transform 的含义

Transform 在 2D/3D 下分别定义为：

• 2D：using Transform = BaseTransform<Rotation2d, Vec2d>，其中 Rotation2d = Eigen::Rotation2D<Real>

• 3D：using Transform = BaseTransform<Rotation3d, Vec3d>，其中 Rotation3d = Eigen::AngleAxis<Real>

Transform 同时包含：

• 旋转（rotation）

• 平移（translation）

并且变换顺序是：先旋转，后平移。

只做平移（最常见）

很多例子写 Transform(translation)，等价于“无旋转 + 平移”。

Vecd translation(1.0, 0.5); // 2D
Transform T(translation);
GeometricShapeBox box(T, halfsize);

同时做旋转 + 平移

当你希望把一个“默认以原点为中心、轴对齐”的几何体旋转后再放到某个位置，可以显式构造旋转并传给 Transform(rotation, translation)。

2D（角度单位为弧度，通常逆时针为正，依赖 Eigen 的约定）：

Real theta = 0.25 * M_PI;
Vecd translation(1.0, 0.5);
Transform T(Rotation(theta), translation); // Rotation 在 2D 即 Rotation2d
GeometricShapeBox rotated_box(T, halfsize);

3D（绕某个单位轴旋转角度 theta）：

Real theta = 0.5 * M_PI;
Vecd axis(0.0, 0.0, 1.0);
Vecd translation(0.0, 1.0, 0.0);

Transform T(Rotation(theta, axis.normalized()), translation); // Rotation 在 3D 即 AngleAxis
GeometricShapeBox rotated_box(T, halfsize);

小提示：并不是所有 Shape 的构造函数都直接暴露 Transform 参数。

• 对于像 GeometricShapeBox 这种“简单几何”，通常直接支持 Transform。

• 对于某些更复杂/外部导入的形状（例如 STL），常见做法是：在建模软件里预先旋转模型，或者使用库中提供的“带变换的包装器/映射器”（如果该形状/几何类型支持）。

案例

tests\2d_examples\test_2d_dambreak

使用创建Box的第一个构造函数：

GeometricShapeBox initial_water_block(Transform(water_block_translation), water_block_halfsize, "WaterBody");

BinaryShapes的布尔操作

BinaryShapes的Binary表示布尔运算，Shapes表示这个类会储存多个subshape，也就是说BinaryShapes用多个subshape做布尔运算。布尔运算目前只提供了add（并集）和sub（差集）两种运算。

用法如下：

op<SubShapeType>(args_to_construct_subshape)

其中op为add或subtract，SubShapeType为欲创建的几何类型，如GeometricShapeBox，args_to_construct_subshape为创建SubShapeType对象需提供的参数。

案例

tests\2d_examples\test_2d_dambreak

大矩形减去小矩形（这里并没有用到ComplexShape的特性，令WallBoundary继承于BinaryShapes足矣）：

class WallBoundary : public ComplexShape // ComplexShape可以换成BinaryShapes
{
  public:
    explicit WallBoundary(const std::string &shape_name) : ComplexShape(shape_name)
    {
        add<GeometricShapeBox>(Transform(outer_wall_translation), outer_wall_halfsize);
        subtract<GeometricShapeBox>(Transform(inner_wall_translation), inner_wall_halfsize);
    }
};
int main() {
    ...
    SolidBody wall_boundary(sph_system, makeShared<WallBoundary>("WallBoundary"));
    ...
}

除了创建一个新的继承于BinaryShapes的类，我们还可以直接创建一个BinaryShapes对象来进行布尔操作：

int main() {
    ...
	auto wall_boundary_ = makeShared<BinaryShapes>("WallBoundary");
    wall_boundary_->add<GeometricShapeBox>(Transform(outer_wall_translation), outer_wall_halfsize);
    wall_boundary_->subtract<GeometricShapeBox>(Transform(inner_wall_translation), inner_wall_halfsize);
    SolidBody wall_boundary(sph_system, wall_boundary_);
    ...
}

用水平集方法进行粒子松弛

对于复杂几何（例如斜坡、圆面），lattice方法难以生成光滑的表面（在particle spacing较大时尤为显著）。此时可以进行粒子松弛，让粒子形成贴体分布。粒子松弛会用到水平集方法进行表面约束（Yu et al, 2023）。

使用时，先用–-relax on选项进行粒子松弛，然后用--reload on读取松弛的结果以进行模拟。

案例

tests\2d_examples\test_2d_flow_around_cylinder

int main() {
    ...
    // Cylinder是一个MultiPolygonShape的派生类，定义在2d_flow_around_cylinder.h
	SolidBody cylinder(sph_system, makeShared<Cylinder>("Cylinder"));
    ...
    cylinder.defineBodyLevelSetShape();
    ...
    (!sph_system.RunParticleRelaxation() && sph_system.ReloadParticles())
        ? cylinder.generateParticles<BaseParticles, Reload>(cylinder.getName())
        : cylinder.generateParticles<BaseParticles, Lattice>();
    ...
    if (sph_system.RunParticleRelaxation()) {
        ... // particle relaxation
        return 0;
    }
    ...
}

我们需要创建一个圆形区域，这确定了粒子松弛的初始形状。然后defineBodyLevelSetShape。relax和reload不能同时进行，如果是relax则按lattice生成粒子（尽管后面还是会随机化）；如果是reload，我们首先应确保bin/reload目录下有name_rld.xml文件，name是body的名字（Cylinder），然后使用generateParticles<BaseParticles, Reload>(name)加载粒子。后面如果要relax的话，会有相应的程序，执行完relax会在bin/reload目录下生成name_rld.xml文件，然后直接return 0。

从任意原始几何映射

这类方法的核心思想是：

• 已经有一个“原始形状”（任意 Shape 派生类，如 MultiPolygonShape、TriangleMeshShapeSTL、自定义的 ComplexShape/BinaryShapes 等）。

• 我们不重新写一个新几何，而是“包装/映射”原始几何的接口（主要是 checkContain() / findClosestPoint()），得到一个新的形状。

InverseShape 和 ExtrudeShape 都定义在 src/shared/geometries/mapping_shape.h。

InverseShape：反转内外

InverseShape<BaseShapeType> 会把“内部/外部”反过来：

• 原形状 BaseShapeType::checkContain(p) 为 true（点在形状内）

• 反转后 InverseShape<BaseShapeType>::checkContain(p) 就为 false（点被认为在形状外）

它的实现非常直接：对 checkContain() 的返回值取反。

典型用途：

• 你想把“形状外部区域”当成约束/边界区域使用（例如用于 NearShapeSurface、StaticConfinement 之类的边界条件）。

案例

tests\2d_examples\test_2d_static_confinement

在二维静态约束边界测试中，用 InverseShape<Triangle> 来把三角形区域的“外部”当成 near-surface 的目标：

NearShapeSurface near_surface_triangle(
    water_block, makeShared<InverseShape<Triangle>>("Triangle"));

ExtrudeShape：按厚度扩张/收缩

ExtrudeShape<BaseShapeType> 会用一个厚度 thickness 对原始形状做“表面偏置”（类似膨胀/腐蚀的效果）：

• thickness > 0：扩张（向外长厚一圈）

• thickness < 0：收缩（向内缩一圈）

它通过原始形状的 findClosestPoint() 计算点到表面的距离，并用 thickness 判断点是否应被视为“在扩张/收缩后的形状内”。

典型用途：

• 给复杂几何生成 level set 前做“圆角化/清理尖角”（对粒子生成和粒子松弛更友好）。

• 用“外扩形状 - 原始形状”构造一个薄壳/缓冲层。

案例 tests\2d_examples\test_2d_particle_generator_single_resolution

下面的写法先把原始多边形 original_logo 外扩一层，再减去原始形状，得到“外扩缓冲层”（可以理解为外边界的一圈壳层）：

MultiPolygon original_logo;
original_logo.addAPolygonFromFile(input_body, ShapeBooleanOps::add);

add<ExtrudeShape<MultiPolygonShape>>(4.0 * resolution_ref, original_logo);
subtract<MultiPolygonShape>(original_logo);

案例

tests\3d_examples\test_3d_particle_relaxation_single_resolution

对 STL 形状同样可以外扩一层再减去原始网格，用于清理尖角/构造壳层：

add<ExtrudeShape<TriangleMeshShapeSTL>>(4.0 * dp_0, full_path_to_file, translation, scaling);
subtract<TriangleMeshShapeSTL>(full_path_to_file, translation, scaling);

二维专用方法

用MultiPolygonShape创建

基于Boost geometry。只能用于二维几何。

MultiPolygonShape具有一个MultiPolygon类型的成员multi_polygon_。基于multi_polygon_，用户可以创建许多几何。并且与BinaryShapes不同，MultiPolygon对象支持四种布尔运算（并集add、差集sub、异或/对称差symm_diff和交集intersect）：

• void addAMultiPolygon(MultiPolygon &multi_polygon, ShapeBooleanOps op);。添加一个MultiPolygon对象，用户需指定布尔运算类型。

• void addABoostMultiPoly(boost_multi_poly &boost_multi_poly, ShapeBooleanOps op);。添加一个Boost MultiPolygon对象，用户需指定布尔运算类型。Boost MultiPolygon对象与Multipolygon对象的区别在于前者是Boost geometry库定义的类型，后者是SPHinXsys自己定义的类型。

• void addAPolygon(const std::vector<Vecd> &points, ShapeBooleanOps op);。添加一个多边形，用户需指定多边形顶点坐标和布尔运算类型。对于闭合的多边形，用户需确保首尾点相同。

• void addABox(Transform transform, const Vecd &halfsize, ShapeBooleanOps op);。添加一个矩形，用户需指定矩形的 transform、半尺寸和布尔运算类型。

  这里的 Transform 不仅能表示平移，也能表示旋转 + 平移。因此除了“把矩形放到某个中心点”，也可以创建旋转矩形。

  // 2D：创建一个绕原点旋转 theta 后，再平移到 translation 的矩形
  Real theta = 15.0 / 180.0 * M_PI;
  Vecd translation(0.5, 0.0);
  Transform T(Rotation(theta), translation);
  multi_polygon_.addABox(T, halfsize, ShapeBooleanOps::add);

• void addACircle(const Vecd &center, Real radius, int resolution, ShapeBooleanOps op);。添加一个圆，用户需指定圆心、半径、分辨率（圆上顶点数）和布尔运算类型。

• void addAPolygonFromFile(std::string file_path_name, ShapeBooleanOps op, Vecd translation = Vecd::Zero(), Real scale_factor = 1.0);。从文件中添加一个多边形，用户需指定文件路径名、布尔运算类型、平移矢量和缩放因子。文件中应有两列数据，分别为x坐标和y坐标，用制表符或空格符分隔，没有表头。

案例

tests\2d_examples\test_2d_channel_flow_fluid_shell

这里用addAPolygon方法添加了一个矩形，矩形有四个顶点，但是为了闭合多边形，在water_block_shape最后添加了一个和第一个点相同的点，一共有五个点。

class WaterBlock : public MultiPolygonShape
{
  public:
    explicit WaterBlock(const std::vector<Vecd> &shape, const std::string &shape_name) : MultiPolygonShape(shape_name)
    {
        multi_polygon_.addAPolygon(shape, ShapeBooleanOps::add);
    }
};
    
void channel_flow_shell(...) {
	...
    /** create a water block shape */
    auto createWaterBlockShape = [&]()
    {
        // geometry
        std::vector<Vecd> water_block_shape;
        water_block_shape.push_back(Vecd(-DL_sponge, 0.0));
        water_block_shape.push_back(Vecd(-DL_sponge, DH));
        water_block_shape.push_back(Vecd(DL, DH));
        water_block_shape.push_back(Vecd(DL, 0.0));
        water_block_shape.push_back(Vecd(-DL_sponge, 0.0));

        return water_block_shape;
    };
	...
    FluidBody water_block(sph_system, makeShared<WaterBlock>(createWaterBlockShape(), "WaterBody"));
}

三维专用方法

用TriangleMeshShape创建

SPHinXsys为我们提供了四种相应的Shape：方块TriangleMeshShapeBrick（类似于GeometricBox）、球体TriangleMeshShapeSphere（类似于GeometricBall）、圆柱TriangleMeshShapeCylinder和基于STL文件创建TriangleMeshShapeSTL。

前三个类基于SimTK库，一个共性参数是resolution，它指定了网格有多精细。注意SimTK的resolution不会影响粒子的间距，只会影响表面有多精细（详见下面的案例），粒子间距是由传入SPHSystem构造函数的resolution_ref控制的。三个类的resolution具体含义各不相同，详见Simbody: SimTK::PolygonalMesh Class Reference。

案例

tests\3d_examples\test_3d_poiseuille_flow_shell

当我们创建TriangleMeshShapeCylinder时，需要传入的参数有轴线方向、底面半径、半长度、SimTK_resolution（如前所述）和圆柱的中心点坐标。

void poiseuille_flow(const Real resolution_ref, const Real resolution_shell, const Real shell_thickness)
{
    ...
    const int SimTK_resolution = 20;
    const Vec3d translation_fluid(0., full_length * 0.5, 0.);
    ...
	auto water_block_shape = makeShared<ComplexShape>("WaterBody");
    water_block_shape->add<TriangleMeshShapeCylinder>(SimTK::UnitVec3(0., 1., 0.), fluid_radius,
                                                      full_length * 0.5, SimTK_resolution,
                                                      translation_fluid);
    ...
}

注意，当前版本的轴线方向矢量应为SimTK::UnitVec3类型，而中心点坐标应为Vecd类型，这倒挺奇怪的，不理解为什么不一致。

补充一点关于“变换”的理解：

• 这里的 translation_fluid 就是圆柱的中心点平移。

• 对圆柱而言，“旋转”通常通过指定轴线方向来表达（例如轴线沿 y 方向：SimTK::UnitVec3(0., 1., 0.)）。

• 对于 TriangleMeshShapeSTL 这类从文件导入的形状，接口往往只提供 translation / scaling，未必提供任意轴旋转；遇到需要旋转 STL 的场景，通常要在导出 STL 前就把模型转好，或者在上层用可用的变换包装方式（如果该版本/该形状支持）。

下面我对比了不同SimTK_resolution的粒子分布。可以看到粒子间距是一样的，大部分粒子坐标也是一样的，不同的是外围（表面）的例子分布。SimTK_resolution==20时，表面更光滑，接近一个圆形；而减小SimTK_resolution至0时，表面很粗糙，接近一个六边形。

At resolution 0 the base is a hexagon with six triangular faces

从SimTK文档我们也可以看到，当resolution是0时，底面被处理为一个六边形。

用ImageShape创建

ImageShape 用一张三维体数据（通常是“距离场 / distance map”）来表示几何体。典型用途：导入复杂几何后，做贴体粒子生成、粒子松弛（level set/distance map约束），以及近表面粒子加密。

因为使用比较少，而且我没咋看明白有啥用，这里不做详细介绍。

案例

tests/3d_examples/test_3d_load_image/load_image.cpp

用空shape自定义粒子坐标

讲完了Shape的所有派生类，再来讲讲空shape。这里空shape指的是不在构造Shape派生类时传递除了名字之外的任何参数。注意，用户不能创建一个Shape基类对象，因为Shape是抽象基类，用户只能创建空的Shape派生类对象。tests中用的比较多的是DefaultShape，它是ComplexShape的别名。在用不到level set时，也可以使用BinaryShapes来创建空shape。

当我们定义一个空shape，意味着后面粒子生成是不依赖于形状的，一般来说也就是我们需要自定义粒子属性。我能想到的有两种方式，一种方式是自己写一个新的ParticleGenerator，另一种是用reload从文件中读取粒子属性。

写一个新的ParticleGenerator

注意要把粒子放在（假想的）网格中心。多用于生成shell的粒子（模板参数为SurfaceParticles和自定义的标签类）。

案例

tests\2d_examples\test_2d_channel_flow_fluid_shell

对于SurfaceParticles来说，在prepareGeometricData()中，用户不仅需要提供粒子坐标，还需要提供体积、法方向、厚度。

/** Particle generator and constraint boundary for shell baffle. */
class WallBoundary;
template <>
class ParticleGenerator<SurfaceParticles, WallBoundary> : public ParticleGenerator<SurfaceParticles>
{
    Real DL_sponge_;
    Real BW_;
    Real resolution_ref_;
    Real wall_thickness_;

  public:
    explicit ParticleGenerator(SPHBody &sph_body, SurfaceParticles &surface_particles,
                               Real resolution_ref, Real wall_thickness)
        : ParticleGenerator<SurfaceParticles>(sph_body, surface_particles),
          DL_sponge_(20 * resolution_ref), BW_(4 * resolution_ref),
          resolution_ref_(resolution_ref), wall_thickness_(wall_thickness) {};
    void prepareGeometricData() override
    {
        auto particle_number_mid_surface = int((DL + DL_sponge_ + 2 * BW_) / resolution_ref_);
        for (int i = 0; i < particle_number_mid_surface; i++)
        {
            Real x = -DL_sponge_ - BW_ + (Real(i) + 0.5) * resolution_ref_;
            // upper wall
            Real y1 = DH + 0.5 * resolution_ref_;
            addPositionAndVolumetricMeasure(Vecd(x, y1), resolution_ref_);
            Vec2d normal_direction_1 = Vec2d(0, 1.0);
            addSurfaceProperties(normal_direction_1, wall_thickness_);
            // lower wall
            Real y2 = -0.5 * resolution_ref_; // lower wall
            addPositionAndVolumetricMeasure(Vecd(x, y2), resolution_ref_);
            Vec2d normal_direction_2 = Vec2d(0, -1.0);
            addSurfaceProperties(normal_direction_2, wall_thickness_);
        }
    }
};
int main() {
    ...
    SolidBody wall_boundary(sph_system, makeShared<DefaultShape>("Wall"));
    ...
    wall_boundary.generateParticles<SurfaceParticles, WallBoundary>(resolution_ref, wall_thickness);
}

从reload文件读取

reload文件是xml格式的文件，一般由粒子松弛过程生成。用户也可以自己写一个reload文件（后面可以考虑更新geoparticle以导出数据到SPHinXsys），它的格式非常简单。文件应命名为name_rld.xml，name是body的名字。文件内容格式为：

<particles>
    <particle Scalar1="0" Vec1="0.1, 0.2" .../>
    <particle Scalar1="1" Vec1="0.1, 0.3" .../>
    ...
</particles>

所有reload文件应当包含以下字段（因为 BaseParticles::registerPositionAndVolumetricMeasureFromReload() 会从 XML 中读取它们）：

字段名	类型	含义
OriginalID	scalar	原始ID
VolumetricMeasure	scalar	体积
Position	vector	坐标

对于SurfaceParticles，还应当包含以下字段：

字段名	类型	含义
NormalDirection	vector	法方向（从流体指向surface）
Thickness	scalar	厚度

程序中如何reload，前面已经讲过，见用水平集方法进行粒子松弛。

AlignedBox

AlignedBox类继承于GeometricBox类，区别在AlignedBox多了一个对齐轴（alignment axis），所有“上/下边界、周期映射、近边界判断”等操作都沿着这个对齐轴来做。

参考文献

Yu, Y., Zhu, Y., Zhang, C., Haidn, O. J. & Hu, X. Level-set based pre-processing techniques for particle methods. Comput. Phys. Commun. 289, 108744 (2023).