二维的三角形网格剖分程序 easymesh.doc

/*% % % o88鞍88o o88鞍88o o88鞍88o 888 888 % 888 888 鞍? 888 888 鞍? 888 888 % 888 888 888 888 888 888 % 8888888? o8888888 ?88888o 888 888 % 888 888 888 888 888 888 % 888 ooo 888 888 ooo 888 ?888888 % ?8oo88? ?888?8 ?8oo88? 888 % ooo 888 % ?88888? % % % 888 % 888 % 888 % 8888888o888o o88鞍88o o88鞍88o 888鞍88o % 888 888 88o 888 888 888 鞍? 888 888 % 888 888 888 888 888 888 888 888 % 888 888 888 8888888? ?88888o 888 888 % 888 888 888 888 888 888 888 % 888 888 888 888 ooo ooo 888 888 888 % 888 888 888 ?8oo88? ?8oo88? 888 888 % % % % Author: Bojan NICENO % % nicenouniv.trieste.it % %*/#define GRAPHICS OFF#include #include #include #include #ifndef max#define max(a,b) (a) (b) ? (a) : (b)#endif#ifndef min#define min(a,b) (a) (b) ? (a) : (b)#endif#ifndef PI#define PI 3endif#define SMALL 1e-30#define GREAT 1e+30#define ON 0 #define OFF -1 /* element is switched off */#define WAIT -2 /* node is waiting (it is a corner node) */#define A 3#define D 4#define W 5#define MAX_NODES 3000/*-+| definitions for the chains |+-*/#define CLOSED 0#define OPEN 1#define INSIDE 2struct ele int i, j, k; int ei, ej, ek; int si, sj, sk; int mark; /* is it off (ON or OFF) */ int state; /* is it (D)one, (A)ctive or (W)aiting */ int material; double xv, yv, xin, yin, R, r, Det; int new_numb; /* used for renumeration */ elemMAX_NODES*2;struct sid int ea, eb; /* left and right element */ int a, b, c, d; /* left, right, start and end point */ int mark; /* is it off, is on the boundary */ double s; int new_numb; /* used for renumeration */ sideMAX_NODES*3;/*sid.ea/eb/a/b/c/d/mark/s/ */struct nod double x, y, F; double sumx, sumy; int Nne; int mark; /* is it off */ int next; /* next node in the boundary chain */ int chain; /* on which chains is the node */ int inserted; int new_numb; /* used for renumeration */ nodeMAX_NODES, pointMAX_NODES/2; struct seg int n0, n1; int N; int chain; int bound; int mark; *segment; struct chai int s0, s1, type; *chain;int Ne, Nn, Ns, Nc; /* number of: elements, nodes, sides */int ugly; /* mora li biti globalna ? */*=*/double area(struct nod *na, struct nod *nb, struct nod *nc) return 0.5 * ( (*nb).x-(*na).x)*(*nc).y-(*na).y) - (*nb).y-(*na).y)*(*nc).x-(*na).x);/*-*/*=*/double dist(struct nod *na, struct nod *nb) return sqrt( (*nb).x-(*na).x)*(*nb).x-(*na).x) + (*nb).y-(*na).y)*(*nb).y-(*na).y) );/*-*/*=*/in_elem(struct nod *n) int e; for(e=0; e= 0.0 & area(n, &nodeeleme.j, &nodeeleme.k) = 0.0 & area(n, &nodeeleme.k, &nodeeleme.i) = 0.0 ) break; return e;/*-in_elem-*/*=*/bowyer(int n, int spac) int e, i, s, swap; struct nod vor; do swap=0; for(s=0; s1 & nodesides.d.bound=OFF & sides.s1 & nodesides.c.bound=OFF & sides.s(nodesides.c.F+nodesides.d.F) ) ) ) */ if(sides.a=n) e=sides.eb; if(e!=OFF)vor.x=eleme.xv; vor.y=eleme.yv; if( dist(&vor, &noden) eleme.R ) swap_side(s); swap=1; else if(sides.b=n) e=sides.ea; if(e!=OFF)vor.x=eleme.xv; vor.y=eleme.yv; if( dist(&vor, &noden) 0) eleme.xv = x = xjk + num/den*(yk-yj); eleme.yv = y = yjk - num/den*(xk-xj); eleme.R = sqrt( (xi-x)*(xi-x) + (yi-y)*(yi-y) ); si=sideeleme.si.s; sj=sideeleme.sj.s; sk=sideeleme.sk.s; O =si+sj+sk; eleme.Det = xi*(yj-yk) - xj*(yi-yk) + xk*(yi-yj); eleme.xin = ( xi*si + xj*sj + xk*sk ) / O; eleme.yin = ( yi*si + yj*sj + yk*sk ) / O; eleme.r = eleme.Det / O;/*-circles-*/*=*/spacing(int e, int n)/*-+| This function calculates the value of the spacing function in | a new node n which is inserted in element e by a linear | approximation from the values of the spacing function in the | elements nodes. |+-*/ double dxji, dxki, dyji, dyki, dx_i, dy_i, det, a, b; dxji = nodeeleme.j.x - nodeeleme.i.x; dyji = nodeeleme.j.y - nodeeleme.i.y; dxki = nodeeleme.k.x - nodeeleme.i.x; dyki = nodeeleme.k.y - nodeeleme.i.y; dx_i = noden.x - nodeeleme.i.x; dy_i = noden.y - nodeeleme.i.y; det = dxji*dyki - dxki*dyji; a = (+ dyki*dx_i - dxki*dy_i)/det; b = (- dyji*dx_i + dxji*dy_i)/det; noden.F = nodeeleme.i.F + a*(nodeeleme.j.F - nodeeleme.i.F) + b*(nodeeleme.k.F - nodeeleme.i.F);/*-spacing-*/*=*/insert_node(double x, double y, int spac, int prev_n, int prev_s_mark, int mark, int next_s_mark, int next_n) int i,j,k,en, n, e,ei,ej,ek, s,si,sj,sk; double sx, sy; Nn+; /* one new node */ nodeNn-1.x = x; nodeNn-1.y = y; nodeNn-1.mark = mark;/* find the element which contains new node */ e = in_elem(&nodeNn-1);/* calculate the spacing function in the new node */ if(spac=ON) spacing(e, Nn-1); i =eleme.i; j =eleme.j; k =eleme.k; ei=eleme.ei; ej=eleme.ej; ek=eleme.ek; si=eleme.si; sj=eleme.sj; sk=eleme.sk; Ne+=2; Ns+=3;/*-+| new elements |+-*/ elemNe-2.i=Nn-1; elemNe-2.j=k; elemNe-2.k=i; elemNe-1.i=Nn-1; elemNe-1.j=i; elemNe-1.k=j; elemNe-2.ei=ej; elemNe-2.ej=Ne-1; elemNe-2.ek=e; elemNe-1.ei=ek; elemNe-1.ej=e; elemNe-1.ek=Ne-2; elemNe-2.si=sj; elemNe-2.sj=Ns-2; elemNe-2.sk=Ns-3; elemNe-1.si=sk; elemNe-1.sj=Ns-1; elemNe-1.sk=Ns-2; /*-+ | new sides |+-*/ sideNs-3.c =k; sideNs-3.d =Nn-1; /* c-d */ sideNs-3.a =j; sideNs-3.b =i; /* a-b */ sideNs-3.ea=e; sideNs-3.eb=Ne-2; sideNs-2.c =i; sideNs-2.d =Nn-1; /* c-d */ sideNs-2.a =k; sideNs-2.b =j; /* a-b */ sideNs-2.ea=Ne-2; sideNs-2.eb=Ne-1; sideNs-1.c =j; sideNs-1.d =Nn-1; /* c-d */ sideNs-1.a =i; sideNs-1.b =k; /* a-b */ sideNs-1.ea=Ne-1; sideNs-1.eb=e; for(s=1; s=3; s+) sx = nodesideNs-s.c.x - nodesideNs-s.d.x; sy = nodesideNs-s.c.y - nodesideNs-s.d.y; sideNs-s.s = sqrt(sx*sx+sy*sy); eleme.i = Nn-1; eleme.ej = Ne-2; eleme.ek = Ne-1; eleme.sj = Ns-3; eleme.sk = Ns-1; if(sidesi.a=i) sidesi.a=Nn-1; sidesi.ea=e; if(sidesi.b=i) sidesi.b=Nn-1; sidesi.eb=e; if(sidesj.a=j) sidesj.a=Nn-1; sidesj.ea=Ne-2; if(sidesj.b=j) sidesj.b=Nn-1; sidesj.eb=Ne-2; if(sidesk.a=k) sidesk.a=Nn-1; sidesk.ea=Ne-1; if(sidesk.b=k) sidesk.b=Nn-1; sidesk.eb=Ne-1; if(ej!=-1) if(elemej.ei=e) elemej.ei=Ne-2; if(elemej.ej=e) elemej.ej=Ne-2; if(elemej.ek=e) elemej.ek=Ne-2; if(ek!=-1) if(elemek.ei=e) elemek.ei=Ne-1; if(elemek.ej=e) elemek.ej=Ne-1; if(elemek.ek=e) elemek.ek=Ne-1;/* Find circumenters for two new elements, and for the one whos segment has changed */ circles(e); circles(Ne-2); circles(Ne-1); bowyer(Nn-1, spac);/*-+| NEW ! Insert boundary conditions for the sides |+-*/ for(s=3; sNs; s+) if(sides.c=prev_n & sides.d=Nn-1) sides.mark=prev_s_mark; if(sides.d=prev_n & sides.c=Nn-1) sides.mark=prev_s_mark; if(sides.c=next_n & sides.d=Nn-1) sides.mark=next_s_mark; if(sides.d=next_n & sides.c=Nn-1) sides.mark=next_s_mark; return e;/*-insert_node-*/*=*/swap_side(int s) int a, b, c, d, ea, eb, eac, ead, ebc, ebd, sad, sac, sbc, sbd; double sx, sy; ea=sides.ea; eb=sides.eb; a=sides.a; b=sides.b; c=sides.c; d=sides.d; if(elemea.ei=eb) ead=elemea.ej; eac=elemea.ek; sad=elemea.sj; sac=elemea.sk; if(elemea.ej=eb) ead=elemea.ek; eac=elemea.ei; sad=elemea.sk; sac=elemea.si; if(elemea.ek=eb) ead=elemea.ei; eac=elemea.ej; sad=elemea.si; sac=elemea.sj; if(elemeb.ei=ea) ebc=elemeb.ej; ebd=elemeb.ek; sbc=elemeb.sj; sbd=elemeb.sk; if(elemeb.ej=ea) ebc=elemeb.ek; ebd=elemeb.ei; sbc=elemeb.sk; sbd=elemeb.si; if(elemeb.ek=ea) ebc=elemeb.ei; ebd=elemeb.ej; sbc=elemeb.si; sbd=elemeb.sj; elemea.i =a; elemea.j =b; elemea.k =d; elemea.ei=ebd; elemea.ej=ead; elemea.ek=eb; elemea.si=sbd; elemea.sj=sad; elemea.sk=s; elemeb.i =a; elemeb.j =c; elemeb.k =b; elemeb.ei=ebc; elemeb.ej=ea; elemeb.ek=eac; elemeb.si=sbc; elemeb.sj=s; elemeb.sk=sac; if(eac!=-1) if(elemeac.ei=ea) elemeac.ei=eb; if(elemeac.ej=ea) elemeac.ej=eb; if(elemeac.ek=ea) elemeac.ek=eb; if(ebd!=-1) if(elemebd.ei=eb) elemebd.ei=ea; if(elemebd.ej=eb) elemebd.ej=ea; if(elemebd.ek=eb) elemebd.ek=ea; if(sidesad.ea=ea) sidesad.a=b; if(sidesad.eb=ea) sidesad.b=b; if(sidesbc.ea=eb) sidesbc.a=a; if(sidesbc.eb=eb) sidesbc.b=a; if(sidesbd.ea=eb) sidesbd.ea=ea; sidesbd.a=a; if(sidesbd.eb=eb) sidesbd.eb=ea; sidesbd.b=a; if(ab) sides.c=a; sides.d=b; sides.a=d; sides.b=c; sides.ea=ea; sides.eb=eb; else sides.c=b; sides.d=a; sides.a=c; sides.b=d; sides.ea=eb; sides.eb=ea; sx = nodesides.c.x - nodesides.d.x; sy = nodesides.c.y - nodesides.d.y; sides.s = sqrt(sx*sx+sy*sy); if(sidesac.ea=ea) sidesac.ea=eb; sidesac.a=b; if(sidesac.eb=ea) sidesac.eb=eb; sidesac.b=b; if(sidesad.ea=ea) sidesad.a=b; if(sidesad.eb=ea) sidesad.b=b; if(sidesbc.ea=eb) sidesbc.a=a; if(sidesbc.eb=eb) sidesbc.b=a; if(sidesbd.ea=eb) sidesbd.ea=ea; sidesbd.a=a; if(sidesbd.eb=eb) sidesbd.eb=ea; sidesbd.b=a; circles(ea); circles(eb);/*-swap_side-*/*=*/void erase() int s, n, e; int a, b, c, d, ea, eb;/*-+| | Negative area check for | elimination of elements | |+-*/ for(e=0; eNe; e+) if( (nodeeleme.i.chain=nodeeleme.j.chain) & (nodeeleme.j.chain=nodeeleme.k.chain) & (chainnodeeleme.i.chain.type=CLOSED) ) a = min( min(eleme.i, eleme.j), eleme.k ); c = max( max(eleme.i, eleme.j), eleme.k ); b = eleme.i+eleme.j+eleme.k - a - c; if(a3) eleme.mark=OFF; else if(area(&nodea, &nodeb, &nodec) 0.0) eleme.mark=OFF; for(e=0; eNe; e+) if(elemeleme.ei.mark=OFF) eleme.ei=OFF; if(elemeleme.ej.mark=OFF) eleme.ej=OFF; if(elemeleme.ek.mark=OFF) eleme.ek=OFF;/*-+| | Elimination of sides | |+-*/ for(s=0; s 3; s+) sides.mark=OFF; for(s=3; sNs; s+) if( (elemsides.ea.mark=OFF) & (elemsides.eb.mark=OFF) ) sides.mark=OFF; for(s=3; sNs; s+) if(sides.mark!=OFF) if(elemsides.ea.mark=OFF) sides.ea=OFF; sides.a=OFF; if(elemsides.eb.mark=OFF) sides.eb=OFF; sides.b=OFF; /*-+| | Elimination of nodes | |+-*/ for(n=0; n 3; n+) noden.mark=OFF;/*-erase-*/*=*/diamond() int ea, eb, eac, ead, ebc, ebd, s; for(s=0; sNs; s+) if(sides.mark!=OFF) ea=sides.ea; eb=sides.eb; if(elemea.ei=eb) ead=elemea.ej; eac=elemea.ek; if(elemea.ej=eb) ead=elemea.ek; eac=elemea.ei; if(elemea.ek=eb) ead=elemea.ei; eac=elemea.ej; if(elemeb.ei=ea) ebc=elemeb.ej; ebd=elemeb.ek; if(elemeb.ej=ea) ebc=elemeb.ek; ebd=elemeb.ei; if(elemeb.ek=ea) ebc=elemeb.ei; ebd=elemeb.ej; if( (eac=OFF | elemeac.state=D) & (ebc=OFF