python三角网格代码_2D三角网格生成程序Distmesh2d

参考资料：

编程及调试遇到的问题和困难：

Delaunay程序主要有以下几点：

1.按给定的信息进行快速查找，需要哈希函数查找模式

2.对三角单元生成进行子女记录，这样插入点在哪个三角单元内部，只需log次操作

3.是否需要考虑插入点在两个相邻三角的交边上。

答：不需要，计算几何书中给出的算法考虑了这点，但编程时却可不考虑，只需对点进行Fuzz操作，即

微小移动，进入到某个三角形内部即可，之后在边翻转操作中会自动处理这个问题，因为这类三角形自身奇异

4.边翻转操作是否有必要模块化：

答：有必要，后期我准备改掉自己的程序，为了加快光滑化网格的速度，现在太慢，Delaunay以OOP的方式编写，并保存信息，由于点的移动不大，我们只需调用子程序对该点翻转，而不需要重复调用Delaunay，大大减少开支

5.随机增量算法，顾名思义，随机加入点是否提高了剖分速度：

答：确实提高了速度，下面的例子会说明。如果不能很好处理点插入边的问题，顺序插入点而非随机插入会更稳健。

关于Distmesh2d的程序：

以论文为依据用Python重写了Distmesh2d程序，网格质量还是不错的，但速度太慢，如何优化，很有必要考虑。遇到的问题如下

1.从Delaunay程序得到的三角信息中，如何提取边

答：我们知道简单的组合边对，会出现大量重复，因为有很多三角形是共用一个边的。

换言之，大部分边，除非边界，都出现过2次，在计算力分量时，需要把重复的去掉。

我的做法是：对每个三角形条每边上2点排序，翻转到对角线以下，利用Python的字典记录出现过的信息，倒序删除。

如果觉得不清楚，请对比下Delauay中的查找

2.如何计算总力

答：作者的论文利用稀疏矩阵存储模式，辅助计算。我的做法是以边进行遍历，计算斥力后，对每个点进行更新。

3.边界拉回操作：作者的程序是1阶精度的，论文里有高阶精度

4.力平衡的自治系统，可以选择更高级的ODE离散模式，论文里说可以隐式求解

下面给出一些Python程序的例子：

以中心为原点，边长为2X2的方形，中间开了0.6半径的圆，运行结果如下：

先看网格图，利用matlab的triplot函数可绘制：

在看下网格的质量，这是内接圆半径与外接圆半径之比的2倍，正三角形为1，退化的则为0

因此1越多越好，那么我们上面的程序的质量如何，在Quality的脚本中，做了计算，并输出到

txt文件，导入matlab绘图请看：

看到这样的网格，还是很欣慰的。有限元分析的误差与网格最小角有关，因此质量近似1的网格，角度近似为

60，这种网格用起来就比较放心对吧

再看时间，这个网格时间代价是高的，由于质量不错，再慢慢优化了，考虑到Delaunay中的随机插入是如何

影响时间，请看下图

26秒和27秒是随机插入，两个32秒是顺序插入，随机插入确实提高了时间，而不是意外如此。这是同一程序，跑了4次

对程序做个小改动，把圆做下偏心，就是这样子

质量不会因此而改变，请看下图，但时间从26秒变到了73秒，我们的点多了一些，时间的增长还是很多的

-----------------------------------------2015-1-15------------------------------

前天，由于网格生成速度太慢，故对其进行了检查，发现Python的数值计算速度太慢。因此用C++对其进行了

重写，如上的网格，C++程序1秒内就给出了答案，而Python要32秒。以图为证：

剩下的问题还有几个：

1.多边形的符号距离函数怎么写，已经有点眉目了；

2.Delauany程序常常因为三点共线退出，需要调整，很迫切；

3.好的初始化网格能大大减少生成时间，针对不同的区域，怎么考究下，发现一个办法是加入某些固定点

对时间的改善明显，比如在圆中心加入一个固定点，等等；

期待后续有更好的发现^ ^

-----------------------------------------------2015-02-16--------------------------------------

任意简单多边形区域的符号距离函数完成了，算法是：

用winding

number判断点在区域内还是在区域外，多边形按旋转顺序依次排序(顺时针或者逆时针)，再计算点到各个相邻线段的投影距离，取绝对值最小者，加上所求符号就是我们的符号距离。

在圆上平均取6个点的多边形就是正六边形，再减去一个圆，得到如下图

如下，C++运行时间2秒

只要细心调整初始化随机删除点的子程序和弹性力平衡距离，就能得到在某些部位加密的效果，

这个可以看作者的论文描述。

C++程序运行时间42秒