24小时热门版块排行榜    

查看: 2901  |  回复: 9

nagami

木虫 (正式写手)

[交流] 2D凸包Delaunay三角剖分程序 已有3人参与

之前用Python写的网格程序,速度实在太慢。没有办法,花了几天时间用C++重写了一次,虽然还有很多瑕疵,不过用用已经可以了,一段实现如下剖分的代码和效果:
资料地址:http://persson.berkeley.edu/distmesh/
程序代码:
#include "nr.h"                          //简易向量矩阵类
#include "distmesh_2d.h"         //参考persson的论文写的基于力平衡的网格程序
#include "delaunay.h"              //随机增量算法Delaunay三角剖分(效果一般)
#include "dfunction.h"            //一般区域的符号距离函数

struct Ftor{
        vector<Point<2>> poly;
        Ftor(vector<Point<2>> &polyy):poly(polyy){}
        Doub fd(const Point<2> &p){
                Doub d1=dcircle(p,0.0,0.0,0.5);                      //圆的符号距离函数
                Doub d2=dpoly(p,poly);                                //多边形的符号距离函数
                return ddiff(d2,d1);                                       //区域做差(挖空圆)
        }
        Doub fh(const Point<2> &p){
                return 0.05+0.3*dcircle(p,0.0,0.0,0.5);           //密度函数,内圆处加密
        }
};
int main(){
        ofstream fp("p.txt";                                              //输出点分布
        ofstream ft("t.txt";                                               //输出三角形节点编号
        Doub pi=3.1415926;                                             
        Box<2> box;                                                       //包围区域初始化剖分点
        box.lo.x[0]=-1.0;box.lo.x[1]=-1.0;
        box.hi.x[0]=1.0;box.hi.x[1]=1.0;

        vector<Point<2>> pfix(6);                               //生成多边形+固定剖分点
        for(int i=0;i<6;i++){                                          //正六变形生成
                Doub theta=pi/3.0;
                pfix.x[0]=cos(theta*i);
                pfix.x[1]=sin(theta*i);
        }
       
        Ftor Fun(pfix);;                                              //距离函数和密度函数初始化
        Distmesh2d<Ftor> mesh(Fun,0.1,box,pfix);  //调用网格程序

      //以下输出,Matlab数组是1偏的,三角形编号+1处理
        fp<<fixed<<setprecision(4);     
        ft<<fixed<<setprecision(6);
        int n=mesh.p.size();
        for(int i=0;i<n;i++){
                fp<<setw(4)<<mesh.p.x[0]<<setw(10)<<mesh.p.x[1]<<endl;
        }
        int k=mesh.tri.size();
        for(int i=0;i<k;i++){
                ft<<setw(10)<<mesh.tri.p[0]+1
                  <<setw(10)<<mesh.tri.p[1]+1
                  <<setw(10)<<mesh.tri.p[2]+1<<endl;
        }
        fp.close();
        ft.close();
        cin.get();
        return 0;
}
2D凸包Delaunay三角剖分程序
密度控制.jpg
2D凸包Delaunay三角剖分程序-1
1.jpg


2D凸包Delaunay三角剖分程序-2
2.jpg



[ Last edited by nagami on 2015-1-16 at 21:30 ]
回复此楼
女靠衣装;男靠金装
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

跑了个随机生成的小圆,缺了个角,我了个擦勒。而且感觉很不靠谱的样子
这可能是什么原因造成的?用matlab的triplot画的,输出脚本应该没问题的。
2D凸包Delaunay三角剖分程序-3
错误了.jpg
女靠衣装;男靠金装
2楼2015-01-08 20:30:02
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

原因找到了:在边翻转操作程序里的指标出现重叠,这种重叠居然造成缺边的现象。
现在看着像合法的三角剖分了
2D凸包Delaunay三角剖分程序-4
P[K8PXQ5`S)EB4LXDOPN6LT.jpg

女靠衣装;男靠金装
3楼2015-01-08 21:57:52
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

人走茶不凉

银虫 (职业作家)


小木虫: 金币+0.5, 给个红包,谢谢回帖
不错,楼主的能力比较强,写的很快!只是好像你随机产生的数据点分布不太均匀。
4楼2015-01-09 00:45:44
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

引用回帖:
4楼: Originally posted by 人走茶不凉 at 2015-01-09 00:45:44
不错,楼主的能力比较强,写的很快!只是好像你随机产生的数据点分布不太均匀。

呵呵,要是产生的网格正确,我也不担心这个。准备在此基础上,做点的平衡分布或者Laplace smoothing,产生高质量的网格。但还有瑕疵,比如插入点在三角形线上就会出现外接圆三点共线的可能问题,会使程序自动推出。
所以再想
1.允许这个瑕疵,因为点的平衡过程会重新移动分布,之后的划分就不会共线。
2.对点进行移动处理,而且还得使其不能跃出最外面的大三角形
3.在原来的基础上添加处理插入点在线上的问题
不知大侠是否研究过此类问题
女靠衣装;男靠金装
5楼2015-01-09 09:11:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

人走茶不凉

银虫 (职业作家)


小木虫: 金币+0.5, 给个红包,谢谢回帖
引用回帖:
5楼: Originally posted by nagami at 2015-01-09 09:11:10
呵呵,要是产生的网格正确,我也不担心这个。准备在此基础上,做点的平衡分布或者Laplace smoothing,产生高质量的网格。但还有瑕疵,比如插入点在三角形线上就会出现外接圆三点共线的可能问题,会使程序自动推出。 ...

我也做过三角剖分的东西,但是也没有太深入研究。
6楼2015-01-10 08:03:54
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

pippi6

铁杆木虫 (著名写手)

工程和科学数值计算咨询


小木虫: 金币+0.5, 给个红包,谢谢回帖
三角网格质量不好,要是能做四边形就好了;三维的能做高质量的hex 就很有用了
7楼2015-01-10 10:52:56
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

2维的凸包Delauany三角剖分程序:Delaunay和基于力平衡网格光滑化程序:Distmesh2d,OOP方式编写,雏形完成了,给自己庆祝下吧
调试了整整两天,各种问题
这里的论文http://persson.berkeley.edu/distmesh/中的代码,和提供的matlab源程序有几处不一样,按论文里进行点的边界拉回,程序失败,因为没注意到这个,一天的时间耗在了这里,逐句对比内部数据的数量级进行排查,来找错,相当耗精力,
正如所见,边界处的毛刺问题没有解决,还有很多优化的地方,如果在秒级的时间内剖分1000左右点的图形,就差不多了。准备花一个礼拜时间慢慢调整
网格预处理结束,就可以好好欣赏有限元理论了
发几张图,留念下
2D凸包Delaunay三角剖分程序-5
边界点的受力扩张移动.jpg


2D凸包Delaunay三角剖分程序-6
程序一角.jpg



2D凸包Delaunay三角剖分程序-7
椭圆剖分.jpg

» 本帖已获得的红花(最新10朵)

女靠衣装;男靠金装
8楼2015-01-11 23:41:26
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

18353248802

捐助贵宾 (小有名气)


小木虫: 金币+0.5, 给个红包,谢谢回帖
能否分享一下MATLAB 的源程序呢?多谢啦
9楼2016-11-26 16:00:15
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

18353248802

捐助贵宾 (小有名气)

送红花一朵
引用回帖:
8楼: Originally posted by nagami at 2015-01-11 23:41:26
2维的凸包Delauany三角剖分程序:Delaunay和基于力平衡网格光滑化程序:Distmesh2d,OOP方式编写,雏形完成了,给自己庆祝下吧
调试了整整两天,各种问题
这里的论文http://persson.berkeley.edu/distme ...

能够分享一下MATLAB源程序呢,多谢啦
10楼2016-11-26 16:01:11
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 nagami 的主题更新
普通表情 高级回复 (可上传附件)
信息提示
请填处理意见