有限元分析中的二维Delaunay三角网格剖分代码实现

有限元分析中的二维Delaunay三角网格剖分代码实现//二维平面点集的Delaunay三角剖分#include "stdafx.h"#include <iostream>#include <GL/glut.h>#include <ctime>#include <cmath>using namespace std;#define point_size 600#define pi 3.1415926struct point{float x,y;};struct triangle{point* Pv[3];float r_of_sqrt;point o_of_tr;};struct d_t_node{triangle Tnode;d_t_node*Pt_l[3];int position_in_dt;int zhuangtai_when_new;};point p1,p2,p3,p4;int n;point p[point_size];int dt_last=0;point p_in_dtriangle1[point_size+1];d_t_node Dtriangle[point_size];point p_in_dtriangle2[point_size+1];d_t_node *queue_t[point_size];point p_in_dtriangle3[point_size+1];int ps_last=0;int queue_t_last=0;point get_spoint_cin(point*p,int n);point get_spoint_rank(point*p,int n);void get_spoint_squre(point*p,int );void get_spoint_cir(point*,int);void read_spoint(point *p,int n);void Display();void init_D_triangle(d_t_node* D_tr);void Pjoin_ps(point*ps,point* p_new,int ps_last);void pjion_D_triangle(d_t_node*Dtriangle,int &dt_last ,point p_new);int findt_clude_p(d_t_node *Dtriangle,int dt_last,point p_new);void sort_to_line(point line_of_tr[point_size][2],int line_of_lin[point_size],int this_fh_when_new[point_size][2], int line_of_tr_last);float hls(float a,float b,float c,float d);point get_o(point p1,point p2,point p3 );float long_of_lines(point point1,point point2);bool p_in_q(d_t_node**queue,point p_new);void put_bug(d_t_node );void put_this(int psize[point_size][2],int last);void read_triangle(d_t_node Dtriangle[point_size] );void put_out_line(point line_of_tr[point_size][2],int line_of_tr_last);int main(int argc, char* argv[]){for(int ii=0;ii<point_size;ii++){for(int jj=0;jj<3;jj++){switch(jj){case 0:Dtriangle[ii].Tnode.Pv[jj]=&p_in_dtriangle1[ii+1];break;case 1:Dtriangle[ii].Tnode.Pv[jj]=&p_in_dtriangle2[ii+1];break;case 2:Dtriangle[ii].Tnode.Pv[jj]=&p_in_dtriangle3[ii+1];break;}}};//点集为圆上的点，输入圆的个数/* int m;cin>>m;n=m*36+1;get_spoint_cir(p,m); *///点集为正方形网格的交点，输入网格的行数/* int m;cin>>m;n=m*m;get_spoint_squre(p,m); *///点集为随机生成的点，输入点的数量cout<<" 请输入点数大小n=";cin>>n;get_spoint_rank(p,n);glutInit(&argc, argv);glutInitDisplayMode(GLUT_RGB | GLUT_SINGLE);glutInitWindowPosition(100, 100);glutInitWindowSize(300, 300);glutCreateWindow("第一个OpenGL程序");glColor3f(0.0,0.0,1.0);init_D_triangle(Dtriangle);for(int i=0;i<n;i++){pjion_D_triangle(Dtriangle,dt_last,p[i]);}glClearColor(1.0,1.0,1.0,1.0);glutDisplayFunc(Display);glutMainLoop();return 0;};void get_spoint_cir(point*p,int nn){float r;p[0].x=0.5;p[0].y=0.5;//int xi,xii;for(int iii=0;iii<nn;iii++ ){//xi=nn-1;xii=xi-iii;r=(float(iii)+1.0)*0.5/float(nn+1);for(int jjj=0;jjj<36;jjj++){p[iii*36+jjj+1].x=float(r*cos(float(jjj)*pi/18.0))+0.5;p[iii*36+jjj+1].y=float(r*sin(float(jjj)*pi/18.0))+0.5;}}};void get_spoint_squre(point*p,int nn){for(int i=0;i<nn;i++)for(int j=0;j<nn;j++)p[i*nn+j].x=float(j+1.0)/float(nn+1),p[i*nn+j].y=float(i+1.0)/float(nn+1);};void put_out_line(point line_of_tr[point_size][2],int line_of_tr_last){for(int i=0;i<line_of_tr_last;i++){cout<<"("<<line_of_tr[i][0].x<<" , "<<line_of_tr[i][0].y<<")-> ";cout<<"("<<line_of_tr[i][1].x<<" , "<<line_of_tr[i][1].y<<") ";}cout<<endl;};void read_triangle(d_t_node Dtriangle[point_size]){for(int i=0;i<dt_last;i++){glBegin(GL_LINE_LOOP);glVertex2f((Dtriangle[i].Tnode.Pv[0]->x-0.5)*2.0,(Dtriangle[i].Tnode.Pv[0]->y-0.5)*2);glVertex2f((Dtriangle[i].Tnode.Pv[1]->x-0.5)*2.0,(Dtriangle[i].Tnode.Pv[1]->y-0.5)*2);glVertex2f((Dtriangle[i].Tnode.Pv[2]->x-0.5)*2.0,(Dtriangle[i].Tnode.Pv[2]->y-0.5)*2);glEnd();};};void read_spoint(point *p,int n){for(int i=0;i<n;i++)glVertex2f(p[i].x,p[i].y);};void Display(){glClear(GL_COLOR_BUFFER_BIT);glPolygonMode(GL_FRONT_AND_BACK,GL_FILL);glPointSize(4.0);read_triangle(Dtriangle);glFlush();}；void put_this(int psize[point_size][2],int last){for(int i=0;i<last;i++)cout<<psize[i][0]<<" ";cout<<endl;for( i=0;i<last;i++)cout<<psize[i][1]<<" ";};void put_bug(d_t_node Dtriangle1){for (int i=0;i<3;i++){cout<<"point : "<<Dtriangle1.Tnode.Pv[i]->x<<" , "<<Dtriangle1.Tnode.Pv[i]->y<<endl;}cout<<"the_o: "<<Dtriangle1.Tnode.o_of_tr.x<<" , "<<Dtriangle1.Tnode.o_of_tr.y<<endl;cout<<"the_r "<<Dtriangle1.Tnode.r_of_sqrt <<endl;cout<<"poistion_in_Dt "<<Dtriangle1.position_in_dt<<endl;cout<<"zhuangtai: "<<Dtriangle1.zhuangtai_when_new<<endl;for (i=0;i<3;i++){if(Dtriangle1.Pt_l[i]!=NULL){cout<<" _ "<<Dtriangle1.Pt_l[i]->position_in_dt;}else cout<<" _NULL ";}cout<<endl;};bool p_in_q(d_t_node*queue_t,point p_new){float p_o,r_t;p_o=long_of_lines(p_new,queue_t->Tnode.o_of_tr);r_t=queue_t->Tnode.r_of_sqrt;if(p_o<=r_t) return 1;else return 0;};float long_of_lines(point point1,point point2){float long_r;long_r=sqrt((point1.x-point2.x)*(point1.x-point2.x)+(point2.y-point1.y)*(point2.y-point1.y)) ;return long_r;};float hls(float a,float b,float c,float d){float e;e=a*d-b*c;return e;};point get_o(point p1,point p2,point p3 ){point o;float a,b,c,d,e;a=(p2.x*p2.x-p3.x*p3.x+p2.y*p2.y-p3.y*p3.y)/2.0;b=p2.y-p3.y;c=(p3.x*p3.x-p1.x*p1.x+p3.y*p3.y-p1.y*p1.y)/2.0;d=p3.y-p1.y;e=hls(p2.x-p3.x,p2.y-p3.y,p3.x-p1.x,p3.y-p1.y);o.x=hls(a,b,c,d)/e;b=(p2.x*p2.x-p3.x*p3.x+p2.y*p2.y-p3.y*p3.y)/2.0;a=p2.x-p3.x;d=(p3.x*p3.x-p1.x*p1.x+p3.y*p3.y-p1.y*p1.y)/2.0;c=p3.x-p1.x;o.y=hls(a,b,c,d)/e;return o;};void sort_to_line(point line_of_tr[point_size][2],int line_of_lin[point_size],int this_fh_when_new[point_size][2], int line_of_tr_last){point change[2];int change1;int change_n1;for(int i=1;i<line_of_tr_last;i++){for(int j=i+1;j<line_of_tr_last;j++){if(line_of_tr[i-1][1].x==line_of_tr[j][0].x&&line_of_tr[i-1][1].y==line_of_tr[j][0].y ) {change[0].x=line_of_tr[j][0].x;change[0].y=line_of_tr[j][0].y;change[1].x=line_of_tr[j][1].x;change[1].y=line_of_tr[j][1].y;line_of_tr[j][0].x=line_of_tr[i][0].x;line_of_tr[j][1].x=line_of_tr[i][1].x;line_of_tr[j][0].y=line_of_tr[i][0].y;line_of_tr[j][1].y=line_of_tr[i][1].y;line_of_tr[i][0].x=change[0].x;line_of_tr[i][0].y=change[0].y;line_of_tr[i][1].x=change[1].x;line_of_tr[i][1].y=change[1].y;change1=this_fh_when_new[i][0];this_fh_when_new[i][0]=this_fh_when_new[j][0];this_fh_ when_new[j][0]=change1;change1=this_fh_when_new[i][1];this_fh_when_new[i][1]=this_fh_when_new[j][1];this_fh_ when_new[j][1]=change1;change_n1=line_of_lin[i];line_of_lin[i]=line_of_lin[j];line_of_lin[j]=change_n1;break;}}}};int findt_clude_p(d_t_node* Dtriangle,int dt_last,point p_new){for(int i=0;i<dt_last;i++){point o_of_t;float r_of_t,r_pando;float a1=Dtriangle[i].Tnode.Pv[0]->x,a2=Dtriangle[i].Tnode.Pv[0]->y;float b1=Dtriangle[i].Tnode.Pv[1]->x,b2=Dtriangle[i].Tnode.Pv[1]->y;float c1=Dtriangle[i].Tnode.Pv[2]->x,c2=Dtriangle[i].Tnode.Pv[2]->y;o_of_t=get_o(*Dtriangle[i].Tnode.Pv[0],*Dtriangle[i].Tnode.Pv[1],*Dtriangle[i].Tnode.Pv[2 ]);r_of_t=sqrt((o_of_t.x-a1)*(o_of_t.x-a1)+(o_of_t.y-a2)*(o_of_t.y-a2));r_pando=sqrt((o_of_t.x-p_new.x)*(o_of_t.x-p_new.x)+(o_of_t.y-p_new.y)*(o_of_t.y-p_new. y));if(r_pando<=r_of_t) return i;}return 0;};void pjion_D_triangle(d_t_node*Dtriangle,int &dt_last,point p_new){queue_t_last=0;d_t_node* t_clude_p=NULL;int here=findt_clude_p(Dtriangle,dt_last, p_new);t_clude_p=&Dtriangle[here];queue_t[queue_t_last]=t_clude_p,queue_t[queue_t_last]->zhuangtai_when_new=1;queue_t_last++;int queue_t_first=0;point line_of_tr[point_size][2];int line_of_tr_last=0;int line_of_lin1[point_size];int k;int this_fh_when_new[point_size][2];while(queue_t_first!=queue_t_last){for(int j=0;j<3;j++){if(queue_t[queue_t_first]->Pt_l[j]!=NULL&&p_in_q(queue_t[queue_t_first]->Pt_l[j],p_new)){if(queue_t[queue_t_first]->Pt_l[j]->zhuangtai_when_new!=1){queue_t[queue_t_last]=queue_t[queue_t_first]->Pt_l[j],queue_t[queue_t_last]->zhuangtai_when_new=1;queue_t_last++;}}else{if(j<2){line_of_tr[line_of_tr_last][0].x=queue_t[queue_t_first]->Tnode.Pv[j]->x;line_of_tr[line_of_tr_last][0].y=queue_t[queue_t_first]->Tnode.Pv[j]->y;line_of_tr[line_of_tr_last][1].x=queue_t[queue_t_first]->Tnode.Pv[j+1]->x;line_of_tr[line_of_tr_last][1].y=queue_t[queue_t_first]->Tnode.Pv[j+1]->y;if(queue_t[queue_t_first]->Pt_l[j]!=NULL)cout<<queue_t[queue_t_first]->Pt_l[j]->position_in_dt<<endl;if(queue_t[queue_t_first]->Pt_l[j]!=NULL){line_of_lin1[line_of_tr_last]=queue_t[queue_t_first]->Pt_l[j]->position_in_dt;for(int l=0;l<3;l++){if(queue_t[queue_t_first]->Pt_l[j]->Pt_l[l]!=NULL&&queue_t[queue_t_first]->Pt_l[j]->Pt_l[ l]->position_in_dt==queue_t[queue_t_first]->position_in_dt)k=l;}this_fh_when_new[line_of_tr_last][0]=queue_t[queue_t_first]->Pt_l[j]->position_in_dt;this_fh_when_new[line_of_tr_last][1]=k;}else{line_of_lin1[line_of_tr_last]=-1;this_fh_when_new[line_of_tr_last][0]=-1;this_fh_when_new[line_of_tr_last][1]=-1;}line_of_tr_last++;}else{line_of_tr[line_of_tr_last][0].x=queue_t[queue_t_first]->Tnode.Pv[j]->x;line_of_tr[line_of_tr_last][0].y=queue_t[queue_t_first]->Tnode.Pv[j]->y;line_of_tr[line_of_tr_last][1].x=queue_t[queue_t_first]->Tnode.Pv[0]->x;line_of_tr[line_of_tr_last][1].y=queue_t[queue_t_first]->Tnode.Pv[0]->y;if(queue_t[queue_t_first]->Pt_l[j]!=NULL){line_of_lin1[line_of_tr_last]=queue_t[queue_t_first]->Pt_l[j]->position_in_dt;for(int l=0;l<3;l++){if(queue_t[queue_t_first]->Pt_l[j]->Pt_l[l]!=NULL&&queue_t[queue_t_first]->Pt_l[j]->Pt_l[l]->p osition_in_dt==queue_t[queue_t_first]->position_in_dt)k=l;}this_fh_when_new[line_of_tr_last][0]=queue_t[queue_t_first]->Pt_l[j]->position_in_dt;this_fh_when_new[line_of_tr_last][1]=k;}else{line_of_lin1[line_of_tr_last]=-1;this_fh_when_new[line_of_tr_last][0]=-1;this_fh_when_new[line_of_tr_last][1]=-1;}line_of_tr_last++;}}}queue_t_first++;}sort_to_line(line_of_tr,line_of_lin1,this_fh_when_new,line_of_tr_last);int last_k,first_k,i;for( i=0;i<queue_t_last;i++){int k=queue_t[i]->position_in_dt;Dtriangle[k].Tnode.Pv[0]->x=p_new.x;Dtriangle[k].Tnode.Pv[0]->y=p_new.y;Dtriangle[k].Tnode.Pv[1]->x=line_of_tr[i][0].x,Dtriangle[k].Tnode.Pv[1]->y=line_of_tr[i][0] .y;Dtriangle[k].Tnode.Pv[2]->x=line_of_tr[i][1].x,Dtriangle[k].Tnode.Pv[2]->y=line_of_tr[i][1] .y;Dtriangle[k].Tnode.o_of_tr=get_o(*Dtriangle[k].Tnode.Pv[0],*Dtriangle[k].Tnode.Pv[1],*Dt riangle[k].Tnode.Pv[2]);Dtriangle[k].Tnode.r_of_sqrt=long_of_lines(Dtriangle[k].Tnode.o_of_tr,*Dtriangle[k].Tnode .Pv[0]);Dtriangle[k].position_in_dt=k;if(i>0){Dtriangle[last_k].Pt_l[2]=&Dtriangle[k];Dtriangle[k].Pt_l[0]=&Dtriangle[last_k];} if(line_of_lin1[i]!=-1){Dtriangle[k].Pt_l[1]=&Dtriangle[line_of_lin1[i]];Dtriangle[this_fh_when_new[i][0]].Pt_l[this_fh_when_new[i][1]]=&Dtriangle[k];}else Dtriangle[k].Pt_l[1]=NULL;Dtriangle[k].zhuangtai_when_new=0;if(i==0)first_k=k;last_k=k;}for(i=0;i<line_of_tr_last-queue_t_last;i++){Dtriangle[i+dt_last].Tnode.Pv[0]->x=p_new.x,Dtriangle[i+dt_last].Tnode.Pv[0]->y=p_new.y;Dtriangle[i+dt_last].Tnode.Pv[1]->x=line_of_tr[i+queue_t_last][0].x,Dtriangle[i+dt_last].Tn ode.Pv[1]->y=line_of_tr[i+queue_t_last][0].y;Dtriangle[i+dt_last].Tnode.Pv[2]->x=line_of_tr[i+queue_t_last][1].x,Dtriangle[i+dt_last].Tn ode.Pv[2]->y=line_of_tr[i+queue_t_last][1].y;Dtriangle[i+dt_last].Tnode.o_of_tr=get_o(*Dtriangle[i+dt_last].Tnode.Pv[0],*Dtriangle[i+dt _last].Tnode.Pv[1],*Dtriangle[i+dt_last].Tnode.Pv[2]);Dtriangle[i+dt_last].Tnode.r_of_sqrt=long_of_lines(Dtriangle[i+dt_last].Tnode.o_of_tr,*Dtri angle[i+dt_last].Tnode.Pv[0]);Dtriangle[i+dt_last].position_in_dt=i+dt_last;if(line_of_lin1[i+queue_t_last]!=-1){Dtriangle[i+dt_last].Pt_l[1]=&Dtriangle[line_of_lin1[i+queue_t_last]];Dtriangle[this_fh_when_new[i+queue_t_last][0]].Pt_l[this_fh_when_new[i+queue_t_last][1] ]=&Dtriangle[i+dt_last];}else Dtriangle[i+dt_last].Pt_l[1]=NULL;if(i>=0){Dtriangle[i+dt_last].Pt_l[0]=&Dtriangle[last_k];Dtriangle[last_k].Pt_l[2]=&Dtriangle[i+dt_last]; }Dtriangle[i+dt_last].zhuangtai_when_new=0;last_k=i+dt_last;}Dtriangle[last_k].Pt_l[2]=&Dtriangle[first_k];Dtriangle[first_k].Pt_l[0]=&Dtriangle[last_k];dt_last=dt_last+line_of_tr_last-queue_t_last;};void Pjoin_ps(point*ps,point* p_new,int ps_last){ps[ps_last]=*p_new;ps_last++;};void init_D_triangle(d_t_node* D_tr){d_t_node *q1=new d_t_node;d_t_node *q2=new d_t_node;D_tr[0].Tnode.Pv[0]->x=0.0; D_tr[0].Tnode.Pv[0]->y=0.0;D_tr[0].Tnode.Pv[1]->x=1.0; D_tr[0].Tnode.Pv[1]->y=0.0;D_tr[0].Tnode.Pv[2]->x=0.0; D_tr[0].Tnode.Pv[2]->y=1.0;D_tr[1].Tnode.Pv[0]->x=1.0; D_tr[1].Tnode.Pv[0]->y=0.0;D_tr[1].Tnode.Pv[1]->x=1.0; D_tr[1].Tnode.Pv[1]->y=1.0;D_tr[1].Tnode.Pv[2]->x=0.0; D_tr[1].Tnode.Pv[2]->y=1.0;D_tr[0].Tnode.o_of_tr=get_o(*D_tr[0].Tnode.Pv[0],*D_tr[0].Tnode.Pv[1],*D_tr[0].Tnode.P v[2]);D_tr[1].Tnode.o_of_tr=get_o(*D_tr[1].Tnode.Pv[0],*D_tr[1].Tnode.Pv[1],*D_tr[1].Tnode.P v[2]);D_tr[0].Tnode.r_of_sqrt=long_of_lines(D_tr[0].Tnode.o_of_tr,*D_tr[0].Tnode.Pv[0]);D_tr[1].Tnode.r_of_sqrt=long_of_lines(D_tr[1].Tnode.o_of_tr,*D_tr[1].Tnode.Pv[0]);D_tr[0].Pt_l[0]=NULL;D_tr[0].Pt_l[1]=&D_tr[1];D_tr[0].Pt_l[2]=NULL;D_tr[1].Pt_l[0]=NULL;D_tr[1].Pt_l[1]=NULL;D_tr[1].Pt_l[2]=&D_tr[0];dt_last=2;D_tr[0].position_in_dt=0;D_tr[1].position_in_dt=1;D_tr[0].zhuangtai_when_new=0, D_tr[1].zhuangtai_when_new=0;};point get_spoint_rank(point*p,int n){point back;back.x=0.0;back.y=0.0;srand(time(0));for(int i=0;i<n;i++){p[i].x=rand()/float(RAND_MAX);if(back.x*back.x<p[i].x*p[i].x) back.x=p[i].x;p[i].y=rand()/float(RAND_MAX);if(back.y*back.y<p[i].x*p[i].y) back.y=p[i].y;cout<<p[i].x<<" "<<p[i].y<<endl;}return back;};point get_spoint_cin(point*p,int n){float M=0,N=0;point back;for(int i=0;i<n;i++){cin>>p[i].x;if(M*M<p[i].x*p[i].x) M=p[i].x;cin>>p[i].y;if(N*N<p[i].y*p[i].y) N=p[i].y;}back.x=M;back.y=N;return back;};。