24小时热门版块排行榜    

查看: 3346  |  回复: 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 的主题更新
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[基金申请] FileCode能看出啥? +8 要乐观耀哥 2026-08-10 23/1150 2026-08-13 00:35 by 要乐观耀哥
[基金申请] 应该是93bebmhtak前后十一个字符比较关键 +22 Lanmanbaby 2026-08-09 36/1800 2026-08-12 22:50 by sdfapple719
[基金申请] 关于Filecode分析方法 +8 majunge000 2026-08-10 9/450 2026-08-12 22:17 by 虫友是什么虫
[基金申请] 好奇怪的filecode +5 布布和一二 2026-08-08 6/300 2026-08-12 16:11 by 云上清扬
[基金申请] filecode +13 documentary 2026-08-10 15/750 2026-08-12 15:56 by 云上清扬
[基金申请] 奇怪,两个人的filecode固定段从头到尾一模一样 +7 布布和一二 2026-08-10 10/500 2026-08-12 09:12 by 布布和一二
[基金申请] 有时候,自然基金真的不能太认真 (我的申报经验) +3 majunge000 2026-08-11 4/200 2026-08-11 20:13 by lch2012
[基金申请] 帮忙看看fileCode +7 wwncly 2026-08-10 13/650 2026-08-11 19:36 by 冰心玉壶晴
[基金申请] 我的国基提前知道中了,可是同事的操作让我实在接受不了,怎么会有这样的人 +8 家与远方 2026-08-10 13/650 2026-08-11 09:49 by 家与远方
[基金申请] 基金中了 +15 laoda193707 2026-08-06 15/750 2026-08-11 00:11 by jiafei2190
[基金申请] 关于代码变化问题,想知道的进来 +17 且听虎啸 2026-08-07 24/1200 2026-08-10 18:35 by zhangduo2008
[基金申请] 国自然结果 +4 Vierhys 2026-08-10 8/400 2026-08-10 15:06 by Vierhys
[基金申请] 据悉今年马上要出结果了 +7 瞬息宇宙 2026-08-10 8/400 2026-08-10 12:42 by Vivilian
[基金申请] 2026国自然放榜时间 +9 布布和一二 2026-08-08 9/450 2026-08-10 11:22 by xxxx2020
[基金申请] 这样的filecode谁见过 +11 布布和一二 2026-08-08 22/1100 2026-08-10 11:10 by wmfsnow
[基金申请] 面上项目filecode邪修 +5 西山十月 2026-08-09 7/350 2026-08-10 07:32 by 仁砚薪传
[基金申请] filecode与中标关系的预测 +5 布布和一二 2026-08-07 5/250 2026-08-09 16:15 by 袁向阳007
[基金申请] 国基金的申报应该改成非等额制,评价高的钱多评价低的钱少,但是增加资助率 +7 a089 2026-08-07 7/350 2026-08-08 18:05 by gltch
[基金申请] 关于filecode +4 布布和一二 2026-08-07 7/350 2026-08-07 22:55 by zhanghaozhu
[基金申请] 大家散了吧,后缀研究没有意义,别浪费时间了,过好目前的每一天,不要焦虑 +5 Tide man 2026-08-06 7/350 2026-08-07 13:11 by 医学老男孩
信息提示
请填处理意见