AutoPacking

1 介绍

实现周期边界，二维任意多边形（凹&凸）和二维圆弧段边界，三维任意多面体（凹&凸的三角片）边界下的填充结构生成。可实现二维圆盘，椭圆，矩形，多边形，圆角多边形等，三维球，椭球，圆柱，球柱，基本柏拉图实体等的填充结构的生成。生成算法包括松弛算法，MonteCarlo算法，DEM算法，RSA算法等。

• 目前开放的部分包括：
  • 松弛算法（球/圆盘）
    • 2d/3d周期边界填充
    • 2d/3d任意多边形/多面体固定边界的填充
    • 3d任意STL三角片边界的填充
    • 2d任意圆弧段边界填充算法

(如需其他算法或者颗粒外形的计算，请联系我们 lsx@pku.edu.cn)

2 使用说明

目前提供三种方式使用该程序，分别是：

• 提供 c++ 的静态链接库，供 c++ 编程使用
• 提供编译好的exe程序，用户编写input.txt输入文件，使用该程序
• 提供基于AutoCAD的插件文件，在CAD中通过相关命令直接使用该程序计算

对于前两种方法，非周期边界还需要准备边界信息文件
对于第三种方法，可直接在通过CAD作图描述边界信息

2.1 边界信息文件说明

• 二维边界
  对于二维边界信息，目前仅仅接受Autocad中pline画出的封闭多边形&圆弧段边界，在cad中画好需要填充的封闭曲线边界后，采用dxfout命令输出边界文件供该程序读取。
  eg: InputFile中的circle.dxf 和 squre10.dxf

• 三维边界
  对于三维边界，目前仅仅接受Autocad中画出的三维实体，采用stlout命令输出的格式文件作为该程序的输入文件。
  eg: InputFile中的box10.stl 和 cylinder.stl 和 sphere8.stl等文件

2.2 c++编程说明

前往AutoPacking网站下载该静态链接库以及相关头文件，在本地解压后，文件目录如下：

AutoPacking.lib
algebra.h
packing.h
packing_monte_carlo.h
packing_relaxation.h
particle.h
particle_circle.h
particle_combined_circle.h
particle_combined_sphere.h
particle_sphere.h
README.md

通过引用所需头文件来实现相关计算
目前已已经开放的的功能：

• 圆盘/球形颗粒
  • 松弛算法
    • 2d/3d周期边界填充算法
    • 2d/3d固定凸边界填充算法
    • 3d固定凹边界填充算法
    • 2d多段线任意边界填充算法
  • 蒙特卡洛算法
    • 2d/3d周期边界填充算法
• 组合圆盘/球
  • 松弛算法
    • 2d/3d周期边界填充算法
  • 蒙特卡洛算法
    • 2d/3d周期边界填充算法

系统要求：windows10 + visual studio 2019
规模限制在10000个颗粒以内，如需更大规模的计算，请联系我们 lsx@pku.edu.cn

2.2.1 软件架构

基类：

particle.h  ///< 颗粒基类
packing.h   ///< 填充算法基类

颗粒的派生类：

particle_circle.h          ///< 圆盘颗粒
particle_combined_circle.h ///< 组合圆盘颗粒
particle_combined_sphere.h ///< 组合球颗粒
particle_sphere.h          ///< 球形颗粒

算法类：

packing_monte_carlo.h   ///< 蒙特卡洛算法
packing_relaxation.h    ///< 填充算法

2.2.2 静态链接库配置教程

windows下配置

1. 在www.packingtech.cn上下载压缩包并本地解压
2. 解决方案配置：Relase
3. 解决方案平台：×64
4. VS->项目->属性-> vc++ 目录->外部包含目录:添加上述目录地址
5. VS->项目->属性-> vc++ 目录->库目录:添加上述目录地址
   目前仅提供vs2019版本和vs2022版本的静态编译库

2.2.3 c++编程说明

算法流程分为以下几个模块
1、定义算法

格式：算法类名<维度> 算法（参数）；
eg:
RelaxationPacking<2> pac(1.5, 1, 1, 0.02, 250);///< 二维松弛算法的定义
RelaxationPacking<3> pac(1.2, 1, 1, 0.02, 250);///< 三维松弛算法
MonteCarloPacking<2> pac(1);///< 二维MC算法的定义
MonteCarloPacking<3> pac(1);///< 三维MC算法的定义
Packing<2> pac;       ///< 定义一个无算法的颗粒容器用于存储2d颗粒

• 其中算法参数为可选项，如无则使用默认参数

2、边界定义
周期边界的定义：

Matrix<2, 2> lattice = 10.0 * Matrix<2, 2>::UNIT();  ///< 定义周期边界矩阵
pac.set_boundary(lattice);                           ///< 将矩阵添加至周期边界

固定边界的定义：

pac.set_boundary("concave_box10.stl");///< 固定边界

• 边界文件的说明同上述文件输入使用该程序

3、设定背景网格类型（此为可选项）

pac.set_cgrid_ratio(2.0);

• 该语句将设定背景网格为Cgrid网格，该网格适用于大小多元填充的计算，相比于NBS背景网格，粒径差别越大，带来的性能提升越大。但是一元情况性能不如NBS背景网格。

4、颗粒生成
定义颗粒：

Particle<2>* circle1 = new CircleParticle(1);定义半径为1的圆盘颗粒
Particle<3>* sphere1 = new SphereParticle(1);定义半径为1的球形颗粒

• <>中的数字表示维度，仅可为2/3，并且圆盘只能是2d颗粒，球只能是3d颗粒。
  将颗粒添加至packing

pac.push(circle1, 500);  ///< 将500个圆盘添加至packing

• 该函数可多次调用以实现多元packing的计算
• 该函数的参数500可省略，则代表仅仅添加一个颗粒
• 该函数在添加多个颗粒时，每个颗粒的坐标为边界内的随机坐标，颗粒方向为随机方向

5、算法初始化

pac.Initialize();  ///< 算法的初始化

• 对于松弛算法，则是初始化颗粒的粒径大小，使之满足初始填充率的设定值
• 对于MC算法，初始化则是初始化颗粒坐标，使颗粒间无重叠，初始化周期边界大小，使之满足初始低密度填充（0.2）

6、开始计算

pac.Run(true);  ///< 开始松弛算法或者MC算法的计算,在屏幕上输出计算过程，此时true可省略
pac.Run(false);///< 开始松弛算法或者MC算法的计算,不在屏幕上输出计算过程

7、结果输出

 pac.ScrPacking("attach_inf",true);                    ///< 输出SCR画图文件
 pac.Output("attach_inf",true);                        ///< 生成xyz格式存储文件

• 其中attach_inf为附加输出的信息，若其后参数为true,则将该信息添加至标准文件名格式之前，若无该参数或者为false,则将该内容添加至输出文件的最后，以文本输出。
• 这两个参数均省略时，则仅仅以标准格式输出文件

2.2.4 2d松弛算法程序示例

///< 二维圆盘松弛程序
#include"particle_circle.h"
#include"packing_relaxation.h"
#pragma comment(lib,"AutoPacking.lib")

int main(void) {
  Matrix<2, 2> lattice = 10.0 * Matrix<2, 2>::UNIT();///< 定义周期边界的大小
  Particle<2>* circle1 = new CircleParticle(1.0);    ///< 定义半径为1.0圆盘
  RelaxationPacking<2> pac(1.5, 1, 1, 0.02, 250);    ///< 二维松弛算法的定
  pac.set_boundary(lattice);        ///< 周期边界
  //pac.set_boundary("squre10.dxf");///< 固定边界
  //pac.set_cgrid_ratio(2.0);       ///< 设置Cgrid计算网格
  pac.push(circle1, 500);           ///< 添加颗粒1至计算程序
  pac.Initialize();                 ///< 松弛计算初始化
  pac.Run();         ///< 松弛计算
  pac.ScrPacking();  ///< 输出SCR画图文件
  pac.Output();      ///< 生成xyz格式存储文件
  return 0;
}

2.3 可执行文件使用说明

前往AutoPacking网站下载该可执行文件，后缀名为 .exe,新建InputFile作为存放输入文件的文件夹，编写input.txt作为参数输入文件，对于非周期边界，还需要边界描述文件，三维情况下为STL文件，二维情况下为DXF文件。
文件目录如下：

InputFile
    input.txt
    boundary.dxf
    boundary.stl
AutoPacking.exe

运行exe文件即可进行计算得到结果

2.3.1 输入文件说明

填充信息输入文件为txt文件

• 第一行为需要计算的packing数量
  对于每一个packing的信息：

• 第一行为packing的维度信息，只允许是2 / 3
• 第二行为填充算法信息：算法名+算法参数

• 松弛算法算法名：relaxation
• 参数：初始填充率 平动松弛系数 转动松弛系数 颗粒缩小系数 松弛周期
• 蒙特卡洛算法名：MonteCarlo
• 参数：压缩速率（每次边界变形每个颗粒的平均移动次数）

• 第三行为边界信息: 边界类型+边界信息

• 边界类型：只允许是 fixed / period
• 固定边界边界信息：文件名，二维填充是 dxf 文件；三维填充是 stl 文件
• 周期边界信息：

• 二维为两个数字，以空格隔开，分别代表周期长方形的长与高
• 三维为三个数字，以空格隔开，分别代表周期长方体的长与宽与高

• 第四行是packing中的组分数量 n ，有多少组分即为多少
• 接下来 n 行是每种组分的信息：包括形状名字，形状参数，该类颗粒数量

• 目前形状名字仅仅包含二维的circle 和三维的sphere
• 对于circle和sphere的形状参数仅仅只有半径R

• 最后一行为分隔符，不具有任何意义，任何连续字符均可

• 示例文件 example_check_algorithm.txt

• 该文件包含13个packing，包括松弛算法固定与周期边界，MC算法固定与周期边界的一元/多元填充算例
• 颗粒类型为圆盘/球
  下面给出一个松弛算法球填充的粒子：

1                               #packing的数量
3                               #维度
relaxation  1.5 1 1 0.005 200   #算法名字和算法参数
fixed box10.stl                 #边界条件
2                               #组分数
sphere 1 200                    #组分1的信息
sphere 1.1 300                  #组分2的信息
****************************    #分隔符

2.4 通过AutoCAD命令使用说明

前往AutoPacking网站下载ObjectARX插件，在AutoCAD中加载该插件后，即可使用
** 张世旋完成后编写该部分 **

2.5 输出文件说明

• 对于每一个packing，均会输出两个文件，分为为基于AutoCAD的scr文件和基于OVITO的xyz格式文件。

2.5.1 文件命名

• 颗粒类型_颗粒数_填充率_填充算法_时间戳.xyz

• 单分散 Sphere(1)[100%]_100_0.641232_MonteCarlo_10011200.xyz

• 形状多分散 Sphere(1)[40%]-Cube(2)[60%]_100_0.765432_Relaxation_10011200.xyz

• 大小多分散 Sphere(1)[40%]-Sphere(1.1)[60%]_100_0.835644_AthermalQuasistatic_10011200.xyz

2.5.2 scr格式说明

• scr文件为AutoCAD中的画图文件，直接用AutoCAD打开即可得到packing图像

2.5.3 xyz格式说明

• 第一行为颗粒数。

• 第二行为周期边界、晶胞中心（固定为原点）、每一列值的含义（用于ovito自动读取）。

• 第三行开始，每一行表示一个颗粒的信息。

  • 其中第一列为颗粒类型，用来区分不同形状、不同大小的颗粒；
  • 后续一列为同形状颗粒类型，用来区分相同形状、不同大小的颗粒；
  • 后续一列为形状名称，取值范围为Shape中的值；
  • 后续三列列为颗粒坐标，分别为x、y、z分量；
  • 后续四列为颗粒方向，以四元数表示，分别为x、y、z、w分量，其中对于圆或球该参数仅作占位；
  • 后续三列为形状轴比，分别为x、y、z分量，其中：
    • 对于圆、球、正方形和立方体，三分量必须相等；
    • 对于多边形、多面体、组合圆（球），该参数仅作占位；
    • 对于圆多边形和球多面体，表示圆角半径；
    • 对于圆（球）柱，x分量表示底面圆半径，z分量表示高度，y分量仅作占位；
  • 最后一列为额外参数，其中第一位表示参数个数，其中：
    • 对于超椭圆（球），表示超椭圆（球）的形状因子；
    • 对于（圆）多边形、（球）多面体，表示局部顶点列表；
    • 对于组合圆（球），表示各圆（球）的半径和局部坐标；

• 整体格式

N
Lattice="R1x R1y R1z R2x R2y R2z R3x R3y R3z" Origin="0.0 0.0 0.0" Properties=species:S:1:local_species:S:1:shape_type:S:1:pos:R:3:orientation:R:4:aspherical_shape:R:3
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
1   0   Ellipse	            Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
2   0   SuperEllipse        Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az	1 p

3   0   Square	            Px	Py	Pz	Qx	Qy	Qz	Qw	L	L	L
4   0   Rectangle           Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
5   0   Polygon	            Px	Py	Pz	Qx	Qy	Qz	Qw	#	#	#	2*n x1 y1 ... xn yn

6   0   CircularPolygon     Px	Py	Pz	Qx	Qy	Qz	Qw	R	R	R	2*n x1 y1 ... xn yn
7   0   CombinedCircle      Px	Py	Pz	Qx	Qy	Qz	Qw	#	#   #   3*n r1 x1 y1 ... rn xn yn

8   0   Sphere    	        Px	Py	Pz	#	#	#	#	R	R	R
9   0   Ellipsoid 	        Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
10  0   SuperEllipsoid      Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az	2 p1 p2

11  0   Cube                Px	Py	Pz	Qx	Qy	Qz	Qw	L	L	L
12  0   Cuboid        	    Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
13  0   Polyhedron    	    Px	Py	Pz	Qx	Qy	Qz	Qw	#	#	#	3*n x1 y1 z1 ... xn yn zn

14  0   SpheroPolyhedron    Px	Py	Pz	Qx	Qy	Qz	Qw	R	R	R	3*n x1 y1 z1 ... xn yn zn
15  0   CombinedSphere      Px	Py	Pz	Qx	Qy	Qz	Qw	#	#	#   4*n r1 x1 y1 z1 ... rn xn yn zn

16  0   Cylinder      	    Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	#	Az
17  0   SpheroCylinder   	Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	#	Az

...

其中，#表示占位符，可置零或任意值，但不可空缺。

• 单分散格式

N
Lattice="R1x R1y R1z R2x R2y R2z R3x R3y R3z" Origin="0.0 0.0 0.0" Properties=species:S:1:local_species:S:1:shape_type:S:1:pos:R:3:orientation:R:4:aspherical_shape:R:3
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
...

• 形状多分散格式

N
Lattice="R1x R1y R1z R2x R2y R2z R3x R3y R3z" Origin="0.0 0.0 0.0" Properties=species:S:1:local_species:S:1:shape_type:S:1:pos:R:3:orientation:R:4:aspherical_shape:R:3
0   0   Circle              Px	Py	Pz	#	#	#	#	R	R	R
1   0   Ellipse	            Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
2   0   Square	            Px	Py	Pz	Qx	Qy	Qz	Qw	L	L	L
3   0   Rectangle           Px	Py	Pz	Qx	Qy	Qz	Qw	Ax	Ay	Az
4   0   CircularPolygon     Px	Py	Pz	Qx	Qy	Qz	Qw	R	R	R	2*n x1 y1 ... xn yn
...

• 大小多分散格式

N
Lattice="R1x R1y R1z R2x R2y R2z R3x R3y R3z" Origin="0.0 0.0 0.0" Properties=species:S:1:local_species:S:1:shape_type:S:1:pos:R:3:orientation:R:4:aspherical_shape:R:3
0   0   Circle              Px	Py	Pz	#	#	#	#	R0	R0	R0
1   1   Circle              Px	Py	Pz	#	#	#	#	R1	R1	R1
2   2   Circle              Px	Py	Pz	#	#	#	#	R2	R2	R2
3   3   Circle              Px	Py	Pz	#	#	#	#	R3	R3	R3
4   4   Circle              Px	Py	Pz	#	#	#	#	R4	R4	R4
...

• 形状大小多分散格式

N
Lattice="R1x R1y R1z R2x R2y R2z R3x R3y R3z" Origin="0.0 0.0 0.0" Properties=species:S:1:local_species:S:1:shape_type:S:1:pos:R:3:orientation:R:4:aspherical_shape:R:3
0   0   Circle              Px	Py	Pz	#	#	#	#	R0	R0	R0
1   1   Circle              Px	Py	Pz	#	#	#	#	R1	R1	R1
2   0   Square	            Px	Py	Pz	Qx	Qy	Qz	Qw	L0	L0	L0
3   1   Square	            Px	Py	Pz	Qx	Qy	Qz	Qw	L1	L1	L1
4   0   CircularPolygon     Px	Py	Pz	Qx	Qy	Qz	Qw	R	R	R	2*n x1 y1 ... xn yn
...

2.5.4 文件注释

除颗粒主要信息外，还需包含

• Periodic: true
• Date: 2021-10-01 12:00:00
• Algorithm: MonteCarlo
• Competation time: 2h34m34s

3 附录

亦可基于该上述提供的静态编译库自行开发相关填充算法，此部分为对于静态库中包含的基类、变量、函数的详细说明

3.1基类说明

particle.hpp和packing.hpp为两个基类
定义新形状需要继承particle.hpp
定义新算法需要继承packing.hpp

• 定义算法的过程中需要考虑对任意形状与组分的通用性

3.2 particle.hpp

目前而言不需要颗粒的材料信息以及运动学信息，因此暂时未定义。
定义新形状必须继承particle.hpp。

3.2.1 数据

 protected:
  Shape shape_type_;         ///< 形状
  Position position_;        ///< 坐标
  Orientation orientation_;  ///< 方向
  int id_;                   ///< 颗粒的id号，与颗粒下标对应

  double volume_;               ///< 体积
  double inner_radius_;         ///< 内切球半径
  double outer_radius_;         ///< 外接球半径
  Vector<D> charactor_length_;  ///<x,y,z方向的颗粒特征长度

基类仅仅包含颗粒的 形状名称，颗粒位置，颗粒方向，颗粒id,内外切圆/球半径，颗粒体积，每个维度上的颗粒特征长度 等基础颗粒信息。

3.2.2 数据操作函数

 public:
  Shape get_shape_type() const ;          ///< 读取形状
  const Position& get_position() const ;  ///< 读取坐标
  const Orientation& get_orientation() const     ///< 读取方向   
  int get_id() const ;            ///<读取颗粒的id
  double get_volume() const ;     ///< 读取体积
  double get_inner_radius() const ;  ///< 读取内切球半径
  double get_outer_radius() const ;  ///< 读取外接球半径
  Vector<D> get_charactor_length() const ; ///<读取颗粒的三个方向的特征长度
 private:
  friend class Packing<D>;
  void set_id(int id) ;  ///<读取颗粒的id
  void set_charactor_length(Vector<D> len) ;    ///<设置颗粒的特征长度
  void set_position(const Position& position);  ///< 设置坐标
  void set_orientation(const Orientation& orientation);  ///< 设置方向

  void adj_position(const Position& delta_position);           ///< 调整坐标
  void adj_orientation(const Orientation& delta_orientation);  ///< 调整方向

说明：
颗粒数据的获取在任何时候均可获取，但是颗粒数据的更改和设置，仅仅可以在Packing类中进行,因为参数改变则需要网格更新同步进行，因此不允许用户调用particle类中的函数设置，仅仅可以调用packing类中的函数更改。

3.2.3 辅助计算函数

virtual void set_particle_scale(double scale) = 0; 

• 将颗粒等比的放大或者缩小，比例因子为scale。
• 由于不同的形状放缩方法与参数可能不同，因此基类中不提供该函数的定义。
• 新定义颗粒必须提供该函数的定义。

virtual std::vector<double> get_shape_parameters()const = 0; 

• 获取除基类中的颗粒参数意外的参数，由于参数个数不定，因此用数组表示。
• 对于圆盘或者球，该参数代表半径。对于椭球，代表每个半长轴的长度。对于更加复杂的形状，该数组每个参数的意义，可由用户自行定义。
• 由于不同的颗粒形状参数不同，因此基类中不提供该函数的定义。
• 新定义颗粒必须提供该函数的定义。

virtual Particle<D>* NewParticle(void) = 0; 

• 生成当前颗粒的一个副本，为其开辟内存空间，用于颗粒体系的生成。
• 返回生成的颗粒的指针。
• 由于不同的颗粒生成方式不同，因此基类中不提供该函数的定义。
• 新定义颗粒必须提供该函数的定义。

virtual bool operator ==(const Particle<DIM>& p) const = 0;

• 定义重载等用于判断两个同名颗粒是否相等
• 比如都是Circle,但是半径不同的两个颗粒应该不相等。
• 新定义颗粒必须提供该函数的定义。

3.2.4 求交函数

• 相交信息的传递

template <int DIM>
struct Collision {
  double depth;         ///< 相交深度
  Vector<DIM> point;    ///< 相交点
  Vector<DIM> normal;   ///< 相交法向
  void Print(void);     ///< 屏幕输出交点信息
};

• 边界数据类型定义如下：

template <int DIM>
struct StlElement {
  Vector<DIM> point[DIM];  ///< DIM==3:三角片的三个顶点,DIM=2:线段的两个顶点
  Vector<DIM> normal;      ///< 边界三角片或者线段的外法向
  void Print(void);        ///< 打印边界信息
  double DistanceToPoint(Vector<DIM> center) const; ///< 计算边界与center 的距离 
};

• 定义颗粒与边界的求交.

  virtual void Overlap(const StlElement<DIM>& stltri,Collision<DIM>& col) const = 0;  

• 颗粒与颗粒求交的中转函数
• 由于是全部使用基类指针表达颗粒，因此该求交函数的作用是通过this指针识别当前颗粒类型后做中转调度使用
• lattice为周期边界下的边界大小，由于我们搜索的都是临近颗粒，当两个颗粒在某一维度上相距过远时（大于两个颗粒的外接球半径之和即可），表明为周期边界，因此在该维度上做周期平移。
• 固定边界的计算该参数可以省略

  virtual void Overlap(const Particle<DIM>* that, Collision<DIM>& col,Matrix<DIM, DIM> lattice = 99999 * Matrix<DIM, DIM>::UNIT())const = 0;  

• 与球/圆形颗粒的求交
• lattice为周期边界下的边界大小，由于我们搜索的都是临近颗粒，当两个颗粒在某一维度上相距过远时（大于两个颗粒的外接球半径之和即可），表明为周期边界，因此在该维度上做周期平移。
• 固定边界的计算该参数可以省略

virtual void Overlap(const SphereParticle* that, Collision<DIM>& col, Matrix<DIM, DIM> lattice = 99999 * Matrix<DIM, DIM>::UNIT()) const ； ///定义颗粒与自身的求交
virtual void Overlap(const CircleParticle* that, Collision<DIM>& col, Matrix<DIM, DIM> lattice = 99999 * Matrix<DIM, DIM>::UNIT()) const ； ///定义颗粒与自身的求交

3.2.5 画图函数

 public:
  virtual void Draw() const;                ///< 自定义绘制图形
  virtual void Print() const;               ///< 屏幕打印输出
  virtual void Input(std::string& line);    ///< 标准格式输入行
  virtual std::string Output() const;       ///< 标准格式输出行
  virtual std::string ScrParticle() const;  ///< scr式输出行

3.3 packing.hpp

• 该类为一个容器，用于存储packing中的颗粒，同时提供自动的背景网格的更新，用于颗粒与颗粒之间的临近检索和颗粒与固定边界之间的临近检索，将全局搜素问题转化为局部搜素问题。
• 该类中包含两套背景网格算法：NBS算法 Cgrid算法
• NBS算法用于一元填充计算，Cgrid算法适用于多元填充计算。

3.3.1 背景网格介绍

• NBS背景网格算法
  如图3.1所示，在2/3维情况下，建立一套固定于坐标系的正方形/立方体背景网格，网格尺寸和颗粒的外接圆/外接球的直径相当。这样可既以保证每一个颗粒最多可以横跨4/8个网格（颗粒质心在网格节点上），又可以保证每个网格内的颗粒尽量最少。每个颗粒所属的网格为其质心所在的网格。
  首先，查找与一个目标颗粒接触的所有颗粒，不需要与系统内所有颗粒进行精确的接触判断，只需要将目标颗粒所在的网格内，以及周围的8/26个网格内的颗粒即可。如图1所示。然后以所有的颗粒作为目标颗粒进行上述操作，就可以找到整个系统内所有的接触。
  

  • 优点：
    • 1）思路简单，便于操作。
    • 2）可以适用于任何维度（1维，2维，3维，甚至更高阶的抽象维数）。
  • 缺点：
    • 1）有大量的重复判断，比如A颗粒和B颗粒的接触，以A作为目标颗粒时判断了一次，以B作为目标颗粒又判断了一次。
      • 求交时可认为减半处理，比如计算颗粒编号为id的颗粒求交时，仅仅计算其与编号大于id的颗粒相交情况，将结果加至两个颗粒即可。
    • 2）对于多相系统，颗粒的形状和大小不一致，尤其是大小差别很大的系统，难以确定一个较好的网格的尺寸。为了保证算法的准确性（只与周围8/26个网格内的颗粒进行判断），网格尺寸必须与最大颗粒的尺寸相当。当颗粒大小差别很大时（如下图），此方法效率很低。
    • 

• Cgrid背景网格介绍
  CGRID算法是基于NBS算法避免重复计算的思路，又考虑颗粒大小不一致的问题的一种算法。具体思路如下：
  CGRID算法中背景网格大小与颗粒尺寸无关，无论多大的背景网格均可保证搜索算法的正确性。在确定颗粒所属的背景网格块时，首先确定颗粒在各个维度上的上下边界，这会在二维时形成一个矩形，在三维形时成一个长方形，与这一包裹颗粒的形状相交的网格都包含了这个颗粒。二维下如下图所示。这样就可以适用于多相系统。
  

  • 优点：
    • 1）有效的减少了重复判断的次数。缩小了判断的范围。
    • 2）可以适用于任何维度（1维，2维，3维，甚至更高阶的抽象维数）。
    • 3）可适用于多相系统，可解决颗粒的形状和大小不一致的问题。
  • 缺点：
    • 1）算法复杂，算法实现比较复杂。
    • 2）网格尺寸可调，不同网格尺寸对效率有很大的影响。最优背景网格尺寸没有确定的取法，需要试算。

3.3.2 packing类中的数据与操作

packing类中包含颗粒数据（对外可见），边界数据（对外可见），填充数据（对外可见）网格数据（对外不可见）。

3.2.2.1 对外可见数据

 private:
  std::vector<Particle<DIM>*> p_particles_;  ///< 颗粒数组
  int num_particles_;                        ///< 记录颗粒的数量
  double particles_volume_;                  ///< 颗粒的总体积
  double max_character_len_;                 ///< 最大的颗粒外切圆半径
  double ave_character_len_;                 ///< 平均颗粒外切圆半径
  std::vector<ShapeNum> packing_information_;  ///< 记录颗粒的形状以及数量

 private:
  Matrix<DIM, DIM> lattice_;               ///< 晶格矩阵
  std::vector<StlElement<DIM>> boundary_;  ///< 边界stl三角片信息
  double grids_volume_;                    ///< 填充空间的总体积
  int is_period_;                   ///< 记录边界是否是周期边界
  double Cgrid_ratio_;    ///< Cgrid网格网格大小比例因子，默认2.0,平均特征长度的2倍
  std::string boundary_file_name_;  ///< 边界信息文件名

  std::string packing_method_name_;  ///<记录填充方法的名字
  bool is_print_;  ///< 是都输出中间计算信息到屏幕

• 以上数据均可以通过函数访问

 public:
  const std::vector<Particle<DIM>*>& get_p_particles() const; ///< 获取颗粒列表
  Particle<DIM>*& operator[](int i);                        ///< 获取单个颗粒
  const std::vector<StlElement<DIM>>& get_boundary() const; ///< 获取边界列表
  const StlElement<DIM>& get_boundary(int id) const;        ///< 获取一条边界
  const std::string get_boundary_file_name() const;         ///< 获取边界文件名
  const double& get_max_character_len() const { return max_character_len_; } ///< 获取packing中的最大特征长度
  const double& get_ave_character_len() const { return ave_character_len_; } ///< 获取packing中的平均特征长度
  const std::vector<ShapeNum>& get_packing_information() const ;  ///< 获取填充信息数据
  const std::string get_packing_method_name() const;  ///< 获取填充方法的名字
  const bool get_is_print() { return is_print_; };    ///< 获得是否输出到屏幕
  const Matrix<DIM, DIM> get_lattice() const ;        ///< 获取模拟域大小，对周期边界即为边界大小
  const bool get_is_period() const;                   ///< 判断是否是周期边界
  const int get_num_particles() const ;      ///< 获取填充中颗粒数
  const double get_grids_volume() const ;    ///< 获取填充中颗粒数量
  const double get_particles_volume() const ;///< 获取颗粒的总体积

• 部分数据可以允许用户通过以下函数设置

 public:
  void set_boundary(std::string filename);           ///< 设置固定边界，该函数只允许调用一次
  void set_boundary(Matrix<DIM, DIM> lattice);       ///< 设置周期边界，该函数只允许调用一次
  void set_lattice(Matrix<DIM, DIM>& lattice);       ///< 设置周期边界的大小
  void set_packing_method_name(std::string method)； ///< 设置方法名
  void set_cgrid_ratio(double ratio)；               ///< 设置cgrid的网格大小，如无该函数，则默认采用NBS网格
  void set_no_print() { is_print_ = false; };        ///< 设置不输出到屏幕

3.2.2.2 对外不可见数据与不可见操作（背景网格相关数据）

 private:
  int grid_num_[DIM];     ///< 每个维度上的网格数量
  double grid_len_[DIM];  ///< 每个维度上的网格长度
  int grid_size_;         ///< 网格总数量
  std::vector<std::array<int, DIM>> grid_pos_;  ///< 储存网格坐标变换（1d->3d）
  bool is_cgrid;  ///< 记录当前背景网格是什么网格，true:Cgrid网格，false:NBS网格

  std::vector<Grid> grids_;  ///< 记录网格中的颗粒以及网格中的边界，以及在NBS中的邻居网格
  std::vector<int> NBS_cell_id_;  ///< 记录颗粒所隶属的网格（仅仅在NBS中使用）
  std::vector<std::vector<int>> Cgrid_cells_id_;  ///< 颗粒所隶属的网格组（仅仅在Cgrid中使用）

其中Grid数据结构如下：

struct Grid {
  std::list<int> particle_indices;  ///< 颗粒数组
  std::vector<int> boundary_indices;  ///<  与该网格中的颗粒可能相交的边界id
  std::vector<int> neighbor_indices;  ///< 邻居列表记录网格中的颗粒以及周围的8/26个网格的编号，仅仅在NBS中使用
};

以上数据仅仅在该类内部可见，也仅仅在该类内部通过函数可操作。

 private:
  void InitialGrid(void);             ///< 网格的初始化
  void AddGridsUpdate(int id);        ///< 在网格中添加颗粒
  void DeleteGridsUpdate(void);       ///< 在网格中删除颗粒
  void ChangeGridsUpdate(int id);     ///< 改变网格中的某个颗粒中的信息
  void ChangeGridsUpdate(void);       ///< 全局更新背景网格
  int pos_to_idx(const std:: array<int, DIM>& grid_pos) const;  ///< 3d to 1d
  std::array<int, DIM> idx_to_pos(int grid_idx) const;         ///< 1d to 3d
  void CheckPeriod(int id);  ///<边界越界检查

• 网格相关处理函数,对外不可见

3.2.2.3 辅助函数与操作

 public:
  ///< 辅助函数
  void Clear();  ///< 清空所有对象
  const double Density() const ;///< 计算填充率
  void ExpandParticles(double scale);  ///< 颗粒放大或者缩小，边界不变
  void ExpandBoundary(double scale);  ///<边界放大或者收缩
  void EliminateOverlap(void);        ///< 将重叠的packing转化为不重叠
  double GetCN(std::vector<int>& cn,double error=0.02);       ///< 计算配位数，输出到数组
  double GetCN(double error = 0.02,std ::string filename = "");  ///< 计算配位数，输出到文件，文件名为空则不输出

  void Translate(const Vector<DIM>& pos, int id);       ///<颗粒平移
  void Rotate(const Orientation& ori, int id);          ///<颗粒旋转
  void SetPosition(const Vector<DIM>& pos, int id);     ///重置颗粒的位置
  void SetOrientation(const Orientation& ori, int id);  ///重置颗粒的方向

  ///< packing 生成相关函数
  void push(Particle<DIM>* particle_template, int size = 1);  ///< 在最后加颗粒，默认一次加一个
  void pop(void);  ///< 删除最后一个颗粒
  void MayOverlapParticle(int id, std::vector<int>& particle_list);  ///<返回可能与颗粒id相交的颗粒id列表
  void MayOverlapBoundary(int id, std::vector<int>& boundary_list);  ///<返回可能相交的颗边界id列表

  ///< 边界相关函数
  Vector<DIM> RandInSTL(void);     ///<在固定边界内部随机生成点
  bool IsInSTL(Vector<DIM>& pos);  ///< 判断点是否在固定边界的内部
  bool InputSTL(std::string filename);  ///<读取固定边界信息
  double CalGridVolume() const;         ///<计算填充空间的体积

  ///<信息输出
  virtual void Print() const;  ///< 屏幕打印信息
  std::string ScrPacking(std::string filename = "",bool isfilename = false);  ///< 打印SCR信息，如果isfilename为true,则将filename输出到文件名，否则输出到文件最后
  void ScrParticle(int id,std::string filename = "");  ///< 打印一个颗粒的scr信息
  std::string Output(std::string filename = "",bool isfilename = false);  ///< 打印packing信息，如果isfilename为true,则将filename输出到文件名，否则输出到文件最后