二维的三角形网格剖分程序 easymeshWord文件下载.docx
《二维的三角形网格剖分程序 easymeshWord文件下载.docx》由会员分享,可在线阅读,更多相关《二维的三角形网格剖分程序 easymeshWord文件下载.docx(62页珍藏版)》请在冰豆网上搜索。
%%%%%%%%%%%%%%%%%%%%%%%%%%niceno@univ.trieste.it%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%*/
#defineGRAPHICSOFF
#include<
math.h>
stdio.h>
stdlib.h>
string.h>
#ifndefmax
#definemax(a,b)(((a)>
(b))?
(a):
(b))
#endif
#ifndefmin
#definemin(a,b)(((a)<
#ifndefPI
#definePI3.14159265359
#defineSMALL1e-30
#defineGREAT1e+30
#defineON0
#defineOFF-1/*elementisswitchedoff*/
#defineWAIT-2/*nodeiswaiting(itisacornernode)*/
#defineA3
#defineD4
#defineW5
#defineMAX_NODES3000
/*-----------------------------+
|definitionsforthechains|
+-----------------------------*/
#defineCLOSED0
#defineOPEN1
#defineINSIDE2
structele
{
inti,j,k;
intei,ej,ek;
intsi,sj,sk;
intmark;
/*isitoff(ONorOFF)*/
intstate;
/*isit(D)one,(A)ctiveor(W)aiting*/
intmaterial;
doublexv,yv,xin,yin,R,r,Det;
intnew_numb;
/*usedforrenumeration*/
}
elem[MAX_NODES*2];
structsid
intea,eb;
/*leftandrightelement*/
inta,b,c,d;
/*left,right,startandendpoint*/
/*isitoff,isontheboundary*/
doubles;
side[MAX_NODES*3];
/*sid.ea/eb/a/b/c/d/mark/s/*/
structnod
doublex,y,F;
doublesumx,sumy;
intNne;
/*isitoff*/
intnext;
/*nextnodeintheboundarychain*/
intchain;
/*onwhichchainsisthenode*/
intinserted;
node[MAX_NODES],point[MAX_NODES/2];
structseg
{intn0,n1;
intN;
intbound;
}
*segment;
structchai
{ints0,s1,type;
*chain;
intNe,Nn,Ns,Nc;
/*numberof:
elements,nodes,sides*/
intugly;
/*moralibitiglobalna?
?
*/
/*=========================================================================*/
doublearea(structnod*na,structnod*nb,structnod*nc)
{
return0.5*(((*nb).x-(*na).x)*((*nc).y-(*na).y)
-((*nb).y-(*na).y)*((*nc).x-(*na).x));
/*-------------------------------------------------------------------------*/
doubledist(structnod*na,structnod*nb)
returnsqrt(((*nb).x-(*na).x)*((*nb).x-(*na).x)
+((*nb).y-(*na).y)*((*nb).y-(*na).y));
in_elem(structnod*n)
inte;
for(e=0;
e<
Ne;
e++)/*Thismustsearchthroughallelements?
if(area(n,&
node[elem[e].i],&
node[elem[e].j])>
=0.0
&
&
area(n,&
node[elem[e].j],&
node[elem[e].k])>
node[elem[e].k],&
node[elem[e].i])>
=0.0)
break;
returne;
/*-in_elem-----------------------------------------------------------------*/
bowyer(intn,intspac)
inte,i,s,swap;
structnodvor;
do
{
swap=0;
for(s=0;
s<
Ns;
s++)
if(side[s].mark==0)
/*if(!
((node[side[s].c].inserted>
1&
node[side[s].d].bound==OFF&
side[s].s<
(node[side[s].c].F+node[side[s].d].F))||
(node[side[s].d].inserted>
node[side[s].c].bound==OFF&
(node[side[s].c].F+node[side[s].d].F))))*/
if(side[s].a==n)
{e=side[s].eb;
if(e!
=OFF)
{vor.x=elem[e].xv;
vor.y=elem[e].yv;
if(dist(&
vor,&
node[n])<
elem[e].R)
{swap_side(s);
swap=1;
}}}
elseif(side[s].b==n)
{e=side[s].ea;
while(swap==1);
/*-bowyer------------------------------------------------------------------*/
circles(inte)
/*---------------------------------------------------+
|Thisfunctioncalculatesradiiofinscribedand|
|circumscribedcircleforagivenelement(inte)|
+---------------------------------------------------*/
doublex,y,xi,yi,xj,yj,xk,yk,xij,yij,xjk,yjk,num,den;
doublesi,sj,sk,O;
xi=node[elem[e].i].x;
yi=node[elem[e].i].y;
xj=node[elem[e].j].x;
yj=node[elem[e].j].y;
xk=node[elem[e].k].x;
yk=node[elem[e].k].y;
xij=0.5*(xi+xj);
yij=0.5*(yi+yj);
xjk=0.5*(xj+xk);
yjk=0.5*(yj+yk);
num=(xij-xjk)*(xj-xi)+(yij-yjk)*(yj-yi);
den=(xj-xi)*(yk-yj)-(xk-xj)*(yj-yi);
if(den>
0)
elem[e].xv=x=xjk+num/den*(yk-yj);
elem[e].yv=y=yjk-num/den*(xk-xj);
elem[e].R=sqrt((xi-x)*(xi-x)+(yi-y)*(yi-y));
si=side[elem[e].si].s;
sj=side[elem[e].sj].s;
sk=side[elem[e].sk].s;
O=si+sj+sk;
elem[e].Det=xi*(yj-yk)-xj*(yi-yk)+xk*(yi-yj);
elem[e].xin=(xi*si+xj*sj+xk*sk)/O;
elem[e].yin=(yi*si+yj*sj+yk*sk)/O;
elem[e].r=elem[e].Det/O;
/*-circles-----------------------------------------------------------------*/
spacing(inte,intn)
/*----------------------------------------------------------------+
|Thisfunctioncalculatesthevalueofthespacingfunctionin|
|anewnode'
n'
whichisinsertedinelement'
e'
byalinear|
|approximationfromthevaluesofthespacingfunctioninthe|
|elementsnodes.|
+----------------------------------------------------------------*/
doubledxji,dxki,dyji,dyki,dx_i,dy_i,det,a,b;
dxji=node[elem[e].j].x-node[elem[e].i].x;
dyji=node[elem[e].j].y-node[elem[e].i].y;
dxki=node[elem[e].k].x-node[elem[e].i].x;
dyki=node[elem[e].k].y-node[elem[e].i].y;
dx_i=node[n].x-node[elem[e].i].x;
dy_i=node[n].y-node[elem[e].i].y;
det=dxji*dyki-dxki*dyji;
a=(+dyki*dx_i-dxki*dy_i)/det;
b=(-dyji*dx_i+dxji*dy_i)/det;
node[n].F=node[elem[e].i].F+
a*(node[elem[e].j].F-node[elem[e].i].F)+
b*(node[elem[e].k].F-node[elem[e].i].F);
/*-spacing-----------------------------------------------------------------*/
insert_node(doublex,doubley,intspac,
intprev_n,intprev_s_mark,intmark,intnext_s_mark,intnext_n)
inti,j,k,en,n,e,ei,ej,ek,s,si,sj,sk;
doublesx,sy;
Nn++;
/*onenewnode*/
node[Nn-1].x=x;
node[Nn-1].y=y;
node[Nn-1].mark=mark;
/*findtheelementwhichcontainsnewnode*/
e=in_elem(&
node[Nn-1]);
/*calculatethespacingfunctioninthenewnode*/
if(spac==ON)
spacing(e,Nn-1);
i=elem[e].i;
j=elem[e].j;
k=elem[e].k;
ei=elem[e].ei;
ej=elem[e].ej;
ek=elem[e].ek;
si=elem[e].si;
sj=elem[e].sj;
sk=elem[e].sk;
Ne+=2;
Ns+=3;
/*---------------+
|newelements|
+---------------*/
elem[Ne-2].i=Nn-1;
elem[Ne-2].j=k;
elem[Ne-2].k=i;
elem[Ne-1].i=Nn-1;
elem[Ne-1].j=i;
elem[Ne-1].k=j;
elem[Ne-2].ei=ej;
elem[Ne-2].ej=Ne-1;
elem[Ne-2].ek=e;
elem[Ne-1].ei=ek;
elem[Ne-1].ej=e;
elem[Ne-1].ek=Ne-2;
elem[Ne-2].si=sj;
elem[Ne-2].sj=Ns-2;
elem[Ne-2].sk=Ns-3;
elem[Ne-1].si=sk;
elem[Ne-1].sj=Ns-1;
elem[Ne-1].sk=Ns-2;
/*------------+
|newsides|
+------------*/
side[Ns-3].c=k;
side[Ns-3].d=Nn-1;
/*c-d*/
side[Ns-3].a=j;
side[Ns-3].b=i;
/*a-b*/
side[Ns-3].ea=e;
side[Ns-3].eb=Ne-2;
side[Ns-2].c=i;
side[Ns-2].d=Nn-1;
side[Ns-2].a=k;
side[Ns-2].b=j;
side[Ns-2].ea=Ne-2;
side[Ns-2].eb=Ne-1;
side[Ns-1].c=j;
side[Ns-1].d=Nn-1;
side[Ns-1].a=i;
side[Ns-1].b=k;
side[Ns-1].ea=Ne-1;
side[Ns-1].eb=e;
for(s=1;
=3;
{sx=node[side[Ns-s].c].x-node[side[Ns-s].d].x;
sy=node[side[Ns-s].c].y-node[side[Ns-s].d].y;
side[Ns-s].s=sqrt(sx*sx+sy*sy);
elem[e].i=Nn-1;
elem[e].ej=Ne-2;
elem[e].ek=Ne-1;
elem[e].sj=Ns-3;
elem[e].sk=Ns-1;
if(side[si].a==i){side[si].a=Nn-1;
side[si].ea=e;
if(side[si].b==i){side[si].b=Nn-1;
side[si].eb=e;
if(side[sj].a==j){side[sj].a=Nn-1;
side[sj].ea=Ne-2;
if(side[sj].b==j){side[sj].b=Nn-1;
side[sj].eb=Ne-2;
if(side[sk].a==k){side[sk].a=Nn-1;
side[sk].ea=Ne-1;
}
if(side[sk].b==k){side[sk].b=Nn-1;
side[sk].eb=Ne-1;
if(ej!
=-1)
{if(elem[ej].ei==e)