24小时热门版块排行榜    

查看: 3354  |  回复: 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

木虫 (正式写手)

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


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



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

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

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

nagami

木虫 (正式写手)

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

nagami

木虫 (正式写手)

原因找到了:在边翻转操作程序里的指标出现重叠,这种重叠居然造成缺边的现象。
现在看着像合法的三角剖分了
2D凸包Delaunay三角剖分程序-7
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的回帖
普通表情 高级回复 (可上传附件)
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 澳大利亚 Murdoch University 全奖博士招生(3个名额)地质化工冶金领域 +13 AI8RaGaPaSCR 2026-08-07 14/700 2026-08-13 17:15 by KC9MWU4rq4LU
[基金申请] 重要来源:本周末出结果 +10 瞬息宇宙 2026-08-12 10/500 2026-08-13 15:46 by likettle
[基金申请] 咱们一起用铁证分析2026国家社科基金中标与否 +3 启萌科技 2026-08-12 6/300 2026-08-13 15:34 by 高飞的梦想
[基金申请] 分享一下我之前已中青C的计划书的filecode +4 布布和一二 2026-08-11 5/250 2026-08-13 12:56 by cratir
[基金申请] 结合人工智能,周易传统文化,filecode打分制来了,3分以上希望很大。 +3 Tide man 2026-08-12 4/200 2026-08-13 08:35 by ZJTJZ
[基金申请] 应该是93bebmhtak前后十一个字符比较关键 +22 Lanmanbaby 2026-08-09 36/1800 2026-08-12 22:50 by sdfapple719
[论文投稿] 职称评审,求友友推荐见刊最快的期刊 +4 工厂打螺丝 2026-08-08 4/200 2026-08-12 09:18 by hansi2025
[基金申请] 是这周出结果还是下周出结果? +3 yuleib84 2026-08-11 3/150 2026-08-11 21:53 by jnhyjjm
[基金申请] 小木虫上这么多卖论文的,真有人买论文么?感觉没必要啊 +8 Tide man 2026-08-10 9/450 2026-08-11 20:39 by beefly
[基金申请] 帮忙看看fileCode +7 wwncly 2026-08-10 13/650 2026-08-11 19:36 by 冰心玉壶晴
[基金申请] 为什么网上很多人说本周 12号出结果 +6 瞬息宇宙 2026-08-10 7/350 2026-08-11 19:25 by Tide man
[基金申请] 什么时候出结果,有咨询渠道??? +3 Tide man 2026-08-11 3/150 2026-08-11 17:54 by kudofaye
[基金申请] 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-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
[基金申请] 化学口download_prp&amp;fileCode的固定段好像这几天一直没变,有变的大神么? +3 Tide man 2026-08-07 4/200 2026-08-07 22:39 by Tide man
[基金申请] 固定端突然变了,今天 +6 archvillain 2026-08-06 10/500 2026-08-07 16:03 by 医学老男孩
[基金申请] filecode变化情况 +6 布布和一二 2026-08-07 22/1100 2026-08-07 14:45 by 且听虎啸
信息提示
请填处理意见