/* pdnmesh - a 2D finite element solver Copyright (C) 2001-2005 Sarod Yatawatta This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA $Id: subdivide.c,v 1.24 2005/02/18 18:27:11 sarod Exp $ */ /* subdivision and refinement of the mesh */ #include #include #include "types.h" /* parameters for modifying the mesh */ MY_DOUBLE g_badness_limit; MY_DOUBLE g_area_floor; MY_DOUBLE g_area_ceil; MY_DOUBLE g_gradient_limit; /* determines line defined by (p1,p2) is * on an edge */ /* returns 0 if no edge, (edge no. +1) if is on an edge */ static int is_on_an_edge(MY_INT p1,MY_INT p2) { MY_INT i; MY_DOUBLE x1,x2,y1,y2; for (i=0; i= 0) && ((My(p1,M)-y2)*(y1-My(p1,M))>= 0) && ((Mx(p2,M)-x2)*(x1-Mx(p2,M))>= 0) && ((My(p2,M)-y2)*(y1-My(p2,M))>= 0) ) { #ifdef DEBUG printf(".. Yes\n"); #endif return(i+1); } } } #ifdef DEBUG printf(".. NO\n"); #endif return(0); } /* a function to calculate the area of a triangle */ static MY_DOUBLE area( triangle *t ) { MY_DOUBLE temp1; temp1 = DETERM(Mx(t->p0,M),My(t->p0,M),Mx(t->p1,M),My(t->p1,M),Mx(t->p2,M),My(t->p2,M)); return(0.5*ABS(temp1)); } /* update neighbour link in a triangle */ /* neighbour old will be replaced by neighbour new */ static void update_neighbours_around(triangle *neighbour,int old,int new) { if ( neighbour ) { if ( neighbour->t0 == old ) { neighbour->t0 = new; } else if ( neighbour->t1 == old ) { neighbour->t1 = new; } else if ( neighbour->t2 == old ) { neighbour->t2 = new; } } } /* a function to determine whether point p4 is inside */ /* the circumcircle of points p1,p2,p3 */ static int inside_circle( int p1, int p2, int p3, int p4) { double temp1,temp2,temp3,temp4,a2,b2; #ifdef DEBUG printf("checking (%d,%d,%d) includes %d\n",p1,p2,p3,p4); #endif temp1 = (Mx(p1,M)*Mx(p1,M)-Mx(p2,M)*Mx(p2,M)+My(p1,M)*My(p1,M)-My(p2,M)*My(p2,M)); temp2 = (Mx(p1,M)*Mx(p1,M)-Mx(p3,M)*Mx(p3,M)+My(p1,M)*My(p1,M)-My(p3,M)*My(p3,M)); temp3 = (Mx(p1,M)-Mx(p2,M))*(My(p1,M)-My(p3,M))-(Mx(p1,M)-Mx(p3,M))*(My(p1,M)-My(p2,M)); /* we know this is not zero */ temp3 = 1/temp3; a2 = (temp1*(My(p1,M)-My(p3,M))-temp2*(My(p1,M)-My(p2,M)))*temp3; b2 = (temp2*(Mx(p1,M)-Mx(p2,M))-temp1*(Mx(p1,M)-Mx(p3,M)))*temp3; temp4 = Mx(p1,M)*Mx(p1,M)-Mx(p4,M)*Mx(p4,M)+My(p1,M)*My(p1,M)-My(p4,M)*My(p4,M); if ( a2*(Mx(p1,M)-Mx(p4,M))+b2*(My(p1,M)-My(p4,M)) < temp4 ) { #ifdef DEBUG printf("points are inside circle\n"); #endif return(1); } return(0); } /* tgn - triangle number to refine */ /* ptn - point number of vertex of tgn facing the edge to be refined */ /* assertion - at each step, the remainder of triangles form a valid * delaunay triangulation */ int delaunay_refine_constrained(int tgn, int ptn, RBT_tree * rbt,DAG_tree *dag) { triangle tq,*t,*t_n; int neighbour,p0,p1,p2,p4; int t1,t2,t3,t4; /* neighbours */ triangle *tp1,*tp2,*tp3,*tp4,*tnew1,*tnew2; DAG_node *current_dag_node,*parent_dag_node; /* fix everything into this generalized form * p2/\ * /| \ * / | \ * t1/ | \ t4 * / | \ * p0 / | \ p4 * \ t | tn / * \ | / * t2 \ | / t3 * \ | / * p1 */ /* get the triangle */ tq.n=tgn; t=RBT_find(&tq, rbt); if ( t && ( t->n == tgn )) { /* triangle exists */ /* first find neighbour facing the edge to be refined */ if ( ptn == t->p0 ) { neighbour=t->t0; p0=t->p0;p1=t->p1;p2=t->p2; t1=t->t1; t2=t->t2; } else if ( ptn == t->p1 ) { neighbour=t->t1; p0=t->p1;p1=t->p2;p2=t->p0; t1=t->t2; t2=t->t0; } else { neighbour=t->t2; p0=t->p2;p1=t->p0;p2=t->p1; t1=t->t0; t2=t->t1; } /* if the edge to be flipped is on a boundary, * we cannot do anything, we quit */ if ( is_on_an_edge(p1,p2) ) { #ifdef DEBUG printf("delaunay_refine_cons: edge on boundary :cannot refine\n"); #endif return(1); } /* now find the neighbouring triangle */ /* it may not exist */ tq.n=neighbour; t_n=RBT_find(&tq, rbt); if ( (!t_n) || (t_n->status!=LIVE) ) {/* no neighbour */ return(0); } /* find the shared edge with neigbour, or rather, the opposing vertex */ if ( tgn == t_n->t0 ) { p4=t_n->p0; t3=t_n->t1; t4=t_n->t2; } else if ( tgn == t_n->t1 ) { p4=t_n->p1; t3=t_n->t2; t4=t_n->t0; } else { p4=t_n->p2; t3=t_n->t0; t4=t_n->t1; } /* if we form a triangle with all three vertices fixed, we cannot refine */ if ((Mval(p0,M)==FX)&&(Mval(p1,M)==FX)&&(Mval(p4,M)==FX)) { #ifdef DEBUG printf("delaunay_refine_cons: fixed points: cannot refine\n"); #endif return(1); } /* now we have fitted everything to the general framework */ if ( inside_circle(p0,p2,p1,p4) ) { /* we flip the edge (p1,p2) */ /* this could only be done if p4 can be seen * by p0, i.e. (p1,p2) is the largest edge of neighbour * but since we assert a valid initial triangulation, * this will hold. */ /* get neighbouring triangles first */ tq.n=t1;tp1=RBT_find(&tq, rbt); tq.n=t2;tp2=RBT_find(&tq, rbt); tq.n=t3;tp3=RBT_find(&tq, rbt); tq.n=t4;tp4=RBT_find(&tq, rbt); /* create 2 new triangles */ if ((tnew1=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tnew2=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } /* we transform the setup into * p2/\ * / \ * / \ * t1/ \ t4 * / new1 \ * p0 / \ p4 * \- - ------ / * \ new2 / * t2 \ / t3 * \ / * p1 */ tnew1->p0=p0;tnew1->p1=p4;tnew1->p2=p2; tnew1->n=ntriangles; tnew2->p0=p0;tnew2->p1=p1;tnew2->p2=p4; tnew2->n=ntriangles+1; tnew1->status=tnew2->status=LIVE; tnew1->boundary=tnew2->boundary=t->boundary; t->status=t_n->status=DEAD; /* neighbours */ tnew1->t0=t4;tnew1->t1=t1;tnew1->t2=tnew2->n; tnew2->t0=t3;tnew2->t1=tnew1->n;tnew2->t2=t2; /* undate the neighbours of the neighbours too */ update_neighbours_around(tp1,tgn,tnew1->n); update_neighbours_around(tp2,tgn,tnew2->n); update_neighbours_around(tp3,neighbour,tnew2->n); update_neighbours_around(tp4,neighbour,tnew1->n); /* insert the 2 new triangles into RB tree */ RBT_insert((void*)tnew1, rbt); parent_dag_node=t->dag_p; current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew1); tnew1->dag_p=current_dag_node; /* link it to the other parent in dag */ DAG_link(dag,t_n->dag_p,current_dag_node); /* second triangle */ RBT_insert((void*)tnew2, rbt); parent_dag_node=t->dag_p; current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew2); tnew2->dag_p=current_dag_node; /* link it to the other parent in dag */ DAG_link(dag,t_n->dag_p,current_dag_node); /* increment triangle count */ ntriangles +=2; /* printf("flipped %d neighbour %d to %d and %d\n",tgn,neighbour,tnew1->n,tnew2->n); */ /* call this function recursively */ delaunay_refine_constrained(tnew1->n,p0,rbt,dag); delaunay_refine_constrained(tnew2->n,p0,rbt,dag); } } else { /* error */ printf("delaunay_refine_cons: cannot find triangle. something is wrong\n"); return(-1); } return(0); } /* decide to split triangle or not */ static int do_we_split_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness) { #ifdef DEBUG printf("do_we_split: badness "MDF">? "MDF", "MDF" > area > "MDF"\n",badness, g_badness_limit,g_area_floor,g_area_ceil); #endif if (( (badness>g_badness_limit) && (tarea>g_area_floor)) || (tarea >g_area_ceil) ) { return(1); } return(0); } /* decide to split triangle or not */ /* also consider gradient */ static int do_we_refine_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness,MY_DOUBLE gradient) { #ifdef DEBUG printf("do_we_split: badness "MDF">? "MDF", "MDF" > area > "MDF" "MDF">"MDF"\n",badness, g_badness_limit,g_area_floor,g_area_ceil,gradient,g_gradient_limit); #endif if ( gradient > g_gradient_limit || ((badness>g_badness_limit) && (tarea>g_area_floor)) || (tarea >g_area_ceil) ) { return(1); } return(0); } /* given triangle t and its neighbour n * determine whether we can split to 4. we can when * they share a boundary edge or the edge * common is the largest of both triangles */ static int can_we_split_to_4(triangle *t, triangle *n) { if(t->boundary != n->boundary ) { return(1); } /* find the largest edges */ /* we know the largest edge of t is one common with n */ /* find the largest edge of n */ if ( n->t0==t->n ) { return((LENGTH(n->p1,n->p2,M)>=LENGTH(n->p0,n->p2,M)) && (LENGTH(n->p1,n->p2,M)>=LENGTH(n->p0,n->p1,M))); } if ( n->t1==t->n ) { return((LENGTH(n->p2,n->p0,M)>=LENGTH(n->p1,n->p2,M)) && (LENGTH(n->p2,n->p0,M)>=LENGTH(n->p1,n->p0,M))); } /* else n->t2==t->n */ return((LENGTH(n->p0,n->p1,M)>=LENGTH(n->p1,n->p2,M)) && (LENGTH(n->p0,n->p1,M)>=LENGTH(n->p0,n->p2,M))); } /* subdivision of a triangle */ int subdivide_this_triangle(triangle *tg) { /* our policy is simple and as follows: * 1) all edges on boundary - insert point at centroid - split to 3 * 2) in all other cases, draw a perpendicular to one of the edges * split to 4 or 2 * the perpendicular is drawn to the longest edge, if other two edges * do not belong to the boundary, else, the free edge is chosen to draw it */ int fixed_points, boundary_edges; point p; DAG_node *current_dag_node,*parent_dag_node; triangle *tg1,*tg2,*tg3,*tg4,*tempt,tq; triangle *t0,*t1,*t2,*t3,*neigh; MY_INT point_number,p0,p1,p2,p3; MY_DOUBLE len0,len1,len2; MY_DOUBLE min_edge,max_edge,triangle_badness; /* how to divide ? */ int how_to_divide; /* 3: divide into 3, inserting point at centroid * 4: divide into 4 (or 2) inserting a normal to longest edge * 0: no division */ /* ignore some warnings about uninitialized variables */ t1=t2=tg3=tg4=0; p2=0; fixed_points=0; if (Mval(tg->p0,M)==FX) { fixed_points++; } if (Mval(tg->p1,M)==FX) { fixed_points++; } if (Mval(tg->p2,M)==FX) { fixed_points++; } boundary_edges=0; if (is_on_an_edge(tg->p0,tg->p1)) { boundary_edges++; } if (is_on_an_edge(tg->p1,tg->p2)) { boundary_edges++; } if (is_on_an_edge(tg->p2,tg->p0)) { boundary_edges++; } /* lengths of edges */ len0=LENGTH(tg->p1,tg->p2,M); len1=LENGTH(tg->p2,tg->p0,M); len2=LENGTH(tg->p1,tg->p0,M); min_edge=(((len0<=len1) && (len0<=len1))? len0 : \ (((len1<=len0) && (len1<=len2))? len1 : len2)); max_edge=(((len0>=len1) && (len0>=len1))? len0 : \ (((len1>=len0) && (len1>=len2))? len1 : len2)); triangle_badness=max_edge/min_edge; #ifdef DEBUG printf("\n\nsubdivide_this: considering triangle %d points(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2); printf("subdivide_this: min edge "MDF", max edge "MDF", badness "MDF" area "MDF"\n",min_edge,max_edge,triangle_badness, area(tg)); printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges); #endif how_to_divide=0; if ( fixed_points == 3 ) { how_to_divide=3; } else if ( do_we_split_this_triangle(area(tg),triangle_badness) ) { /* find neighbour */ if ( (len2 >= len1) && (len2>=len0) ) { tq.n=tg->t2; neigh=RBT_find(&tq,&rt); } else if ((len0 >= len1) && (len0>=len2) ) { tq.n=tg->t0; neigh=RBT_find(&tq,&rt); } else { tq.n=tg->t1; neigh=RBT_find(&tq,&rt); } if ( (!neigh) || (neigh->status!=LIVE ) ) { /* no neighbour exists, divide into 4 or 2 */ how_to_divide=4; } else if ( neigh &&can_we_split_to_4(tg, neigh) ) { /* if the largest edge is a boundary edge * we know it cannot be flipped by delaunay refinement. * so we split to 4 * also if the larges edge of neighbour is the largest of triangle, * we split to 4 */ how_to_divide=4; } else { how_to_divide=3; } /* if the neighbours largest edge is also our largest edge, * we can split to 4 */ } if ( how_to_divide==3 ) { #ifdef DEBUG printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges); #endif /* split to 3 by inserting point at centroid */ /* p0 /\ p0 /|\ * / \ / | \ * / \ to /1 | \ t1 * / \ t2 / / \3 \ * / \ / / pn\ \ * / \ / / \ \ * / p1 \ p2 // 2 \\ * -------------- --------------- * p1 t0 p2 */ p.x=(Mx(tg->p0,M)+Mx(tg->p1,M)+Mx(tg->p2,M))*ONE_OVER_THREE; p.y=(My(tg->p0,M)+My(tg->p1,M)+My(tg->p2,M))*ONE_OVER_THREE; if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) { fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } p.z[0]=(Mz(tg->p0,M)+Mz(tg->p1,M)+Mz(tg->p2,M))*ONE_OVER_THREE; p.val=VAR; point_number=BIT_insert(&M,p); if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles; tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1; tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2; tg1->boundary=tg2->boundary=tg3->boundary=tg->boundary; /* arrange neighbours */ tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2; tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0; tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1; /* new triangles are alive */ tg1->status=tg2->status=tg3->status=LIVE; /* old triangle is dead */ tg->status=DEAD; #ifdef DEBUG printf("subdivide_this: split to %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d)\n",tg1->n,tg1->p0,tg1->p1,tg1->p2,tg2->n,tg2->p0,tg2->p1,tg2->p2,tg3->n,tg3->p0,tg3->p1,tg3->p2); #endif /* update links from neighbours too */ tq.n=tg->t1;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg3->n); #endif update_neighbours_around(tempt,tg->n,tg3->n); tq.n=tg->t2;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg1->n); #endif update_neighbours_around(tempt,tg->n,tg1->n); tq.n=tg->t0;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg2->n); #endif update_neighbours_around(tempt,tg->n,tg2->n); parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); tg2->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); tg3->dag_p=current_dag_node; RBT_insert((void*)tg1, &rt); RBT_insert((void*)tg2, &rt); RBT_insert((void*)tg3, &rt); ntriangles+=3; /* since all edges are on boundary, we cannot refine them */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); delaunay_refine_constrained(tg3->n,point_number,&rt,&dt); return(1); } else if ( how_to_divide==4 ) { #ifdef DEBUG printf("subdivide_this: fixed_points=%d boundaries=%d badness="MDF" are="MDF"\n",fixed_points,boundary_edges,triangle_badness,area(tg)); #endif /* split the largest edge by drawing a perpendicular */ /* general form */ /* p0 p0 * /\ /\ * / \ t3 / |\ * t0 / \ / | \ * / t \ / 1| 2\ * p1 /____p___\ p3 p1/____|___\ p3 * \ / \ | 4 / * \ nei / \ 3| / * \ / \ | / * t1 \ / t2 \ |/ * \/ \/ * p2 p2 */ /* p1-p3 is the largest edge of triangle t */ /* if p1-p3 is a boundary, neighbour is hidden and ignored */ if ( (len2 >= len1) && (len2>=len0) ) { /* edge p0-p1 */ p0=tg->p2; p1=tg->p0; p3=tg->p1; tq.n=tg->t2; neigh=RBT_find(&tq,&rt); tq.n=tg->t1; t0=RBT_find(&tq,&rt); tq.n=tg->t0; t3=RBT_find(&tq,&rt); } else if ((len0 >= len1) && (len0>=len2) ) { /* edge p1-p2 */ p0=tg->p0; p1=tg->p1; p3=tg->p2; tq.n=tg->t0; neigh=RBT_find(&tq,&rt); tq.n=tg->t2; t0=RBT_find(&tq,&rt); tq.n=tg->t1; t3=RBT_find(&tq,&rt); } else { /* edge p2-p0 */ p0=tg->p1; p1=tg->p2; p3=tg->p0; tq.n=tg->t1; neigh=RBT_find(&tq,&rt); tq.n=tg->t0; t0=RBT_find(&tq,&rt); tq.n=tg->t2; t3=RBT_find(&tq,&rt); } /* determine foot of normal to longest edge */ /* note the denominator will be zero only if the triangle is * a right angled one. but then the longest edge will be the diagonal * hence in practice it will never be zero unless points coincide */ /* len0=-((Mx(p1)-Mx(p3))*(Mx(p3)-Mx(p0))+(My(p1)-My(p3))*(My(p3)-My(p0))) /((Mx(p1)-Mx(p0))*(Mx(p1)-Mx(p3))+(My(p1)-My(p0))*(My(p1)-My(p3))); p.x=(len0*Mx(p1)+Mx(p3))/(1+len0); p.y=(len0*My(p1)+My(p3))/(1+len0); */ p.x=(Mx(p1,M)+Mx(p3,M))*0.5; p.y=(My(p1,M)+My(p3,M))*0.5; if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) { fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } /* p.z[0]=(len0*Mz(p1)+Mz(p3))/(1+len0); */ p.z[0]=(Mz(p1,M)+Mz(p3,M))*0.5; point_number=is_on_an_edge(p1,p3); /* point_number= boundary_edge+1, if (p1,p3) lies on an edge */ p.val=(point_number && edge_array[point_number-1].type==DR?FX:VAR); #ifdef DEBUG if (p.val==FX) { printf("subdivide_this: point fixed\n"); } else { printf("subdivide_this: point variable\n"); } #endif point_number=BIT_insert(&M,p); /* check for neighbour */ if (neigh && neigh->status==LIVE) { if (neigh->t0==tg->n) { p2=neigh->p0; tq.n=neigh->t1; t1=RBT_find(&tq,&rt); tq.n=neigh->t2; t2=RBT_find(&tq,&rt); } else if (neigh->t1 == tg->n) { p2=neigh->p1; tq.n=neigh->t2; t1=RBT_find(&tq,&rt); tq.n=neigh->t0; t2=RBT_find(&tq,&rt); } else { p2=neigh->p2; tq.n=neigh->t0; t1=RBT_find(&tq,&rt); tq.n=neigh->t1; t2=RBT_find(&tq,&rt); } #ifdef DEBUG printf("subdivide_this: neighbour %d with neighbours (%d(%d),%d(%d),%d(%d))\n",neigh->n,neigh->t0,neigh->p0,neigh->t1,neigh->p1,neigh->t2,neigh->p2); printf("sundivide_this: generalized to p0=%d,p1=%d,p2=%d,p3=%d\n",p0,p1,p2,p3); #endif } /*split to 4 or 2 */ if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ( neigh && neigh->status==LIVE) { if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } } /* update triangle data */ tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles; tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1; tg1->boundary=tg2->boundary=tg->boundary; /* alive, dead triangles */ tg1->status=tg2->status=LIVE; tg->status=DEAD; if (neigh && neigh->status==LIVE) { tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2; tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3; tg3->status=tg4->status=LIVE; tg3->boundary=tg4->boundary=neigh->boundary; /* at the last moment this */ neigh->status=DEAD; tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1); tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n; tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n; tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n; update_neighbours_around(t0,tg->n,tg1->n); update_neighbours_around(t3,tg->n,tg2->n); update_neighbours_around(t1,neigh->n,tg3->n); update_neighbours_around(t2,neigh->n,tg4->n); /* insert triangles into RBT */ RBT_insert((void*)tg1,&rt); RBT_insert((void*)tg2,&rt); RBT_insert((void*)tg3,&rt); RBT_insert((void*)tg4,&rt); ntriangles+=4; /* update DAG */ parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2); tg2->dag_p=current_dag_node; parent_dag_node=neigh->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg3); tg3->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg4); tg4->dag_p=current_dag_node; #ifdef DEBUG printf("split to 4 triangles %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d)\n",tg1->n,tg1->p0,tg1->p1,tg1->p2,tg2->n,tg2->p0,tg2->p1,tg2->p2,tg3->n,tg3->p0,tg3->p1,tg3->p2,tg4->n,tg4->p0,tg4->p1,tg4->p2); #endif /* refine the triangulation */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); delaunay_refine_constrained(tg3->n,point_number,&rt,&dt); delaunay_refine_constrained(tg4->n,point_number,&rt,&dt); } else { tg1->t0=(neigh?neigh->n:-1);tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1); tg2->t0=(neigh?neigh->n:-1);tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n; update_neighbours_around(t0,tg->n,tg1->n); update_neighbours_around(t3,tg->n,tg2->n); /* insert triangles into RBT */ RBT_insert((void*)tg1,&rt); RBT_insert((void*)tg2,&rt); ntriangles+=2; /* update DAG */ parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2); tg2->dag_p=current_dag_node; /* refine the triangulation */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); } return(1); } return(0); } /* subdivides all live triangles */ int subdivide_all_triangles(void) { triangle *rec; DAG_traverse_list_reset(&dt); rec=DAG_traverse_prune_list(&dt); while(rec){ if (rec &&(rec->status==LIVE)) { subdivide_this_triangle(rec); } rec=DAG_traverse_prune_list(&dt); } return(0); } /* refinement of a triangle */ /* similar to sundivide_this_triangle * except we consired gradient of potentials */ /* if override==1, we split it no matter what */ int refine_this_triangle(triangle *tg, int override,MY_DOUBLE priority) { /* our policy is simple and as follows: * 1) all edges on boundary - insert point at centroid - split to 3 * 2) in all other cases, draw a perpendicular to one of the edges * split to 4 or 2 * the perpendicular is drawn to the longest edge, if other two edges * do not belong to the boundary, else, the free edge is chosen to draw it */ int fixed_points, boundary_edges; point p; DAG_node *current_dag_node,*parent_dag_node; triangle *tg1,*tg2,*tg3,*tg4,*tempt,tq; triangle *t0,*t1,*t2,*t3,*neigh; MY_INT point_number,p0,p1,p2,p3; MY_DOUBLE len0,len1,len2; MY_DOUBLE min_edge,max_edge,triangle_badness,triangle_gradient; /* how to divide ? */ int how_to_divide; /* 3: divide into 3, inserting point at centroid * 4: divide into 4 (or 2) inserting a normal to longest edge * 0: no division */ /* ignore some warnings about uninitialized variables */ t1=t2=tg3=tg4=0; p2=0; fixed_points=0; if (Mval(tg->p0,M)==FX) { fixed_points++; } if (Mval(tg->p1,M)==FX) { fixed_points++; } if (Mval(tg->p2,M)==FX) { fixed_points++; } boundary_edges=0; if (is_on_an_edge(tg->p0,tg->p1)) { boundary_edges++; } if (is_on_an_edge(tg->p1,tg->p2)) { boundary_edges++; } if (is_on_an_edge(tg->p2,tg->p0)) { boundary_edges++; } /* lengths of edges */ len0=LENGTH(tg->p1,tg->p2,M); len1=LENGTH(tg->p2,tg->p0,M); len2=LENGTH(tg->p1,tg->p0,M); min_edge=(((len0<=len1) && (len0<=len1))? len0 : \ (((len1<=len0) && (len1<=len2))? len1 : len2)); max_edge=(((len0>=len1) && (len0>=len1))? len0 : \ (((len1>=len0) && (len1>=len2))? len1 : len2)); triangle_badness=max_edge/min_edge; triangle_badness*=priority; /* mean triangle gradient */ triangle_gradient=ABS(Mz(tg->p0,M)-Mz(tg->p1,M)) + ABS(Mz(tg->p1,M)-Mz(tg->p2,M)) + ABS(Mz(tg->p2,M)-Mz(tg->p0,M)); /* multiply with priority value * priority will change with the number of iterations */ /*triangle_gradient *=priority*0.3333; */ #ifdef DEBUG printf("\n\nrefine_this: considering triangle %d points(%d,%d,%d) priority "MDF"\n",tg->n,tg->p0,tg->p1,tg->p2,priority); printf("refine_this: min edge "MDF", max edge "MDF", badness "MDF" area "MDF" gradient "MDF"\n",min_edge,max_edge,triangle_badness, area(tg),triangle_gradient); printf("refine_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges); #endif how_to_divide=0; if ( fixed_points == 3 ) { how_to_divide=3; } else if ( override || do_we_refine_this_triangle(area(tg),triangle_badness,triangle_gradient) ) { /* find neighbour */ if ( (len2 >= len1) && (len2>=len0) ) { tq.n=tg->t2; neigh=RBT_find(&tq,&rt); } else if ((len0 >= len1) && (len0>=len2) ) { tq.n=tg->t0; neigh=RBT_find(&tq,&rt); } else { tq.n=tg->t1; neigh=RBT_find(&tq,&rt); } if ( (!neigh) || (neigh->status!=LIVE ) ) { /* no neighbour exists, divide into 4 or 2 */ how_to_divide=4; } else if ( neigh &&can_we_split_to_4(tg, neigh) ) { /* if the largest edge is a boundary edge * we know it cannot be flipped by delaunay refinement. * so we split to 4 */ how_to_divide=4; } else { how_to_divide=3; } /* if the neighbours largest edge is also our largest edge, * we can split to 4 */ } if ( how_to_divide==3 ) { #ifdef DEBUG printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges); #endif /* split to 3 by inserting point at centroid */ /* p0 /\ p0 /|\ * / \ / | \ * / \ to /1 | \ t1 * / \ t2 / / \3 \ * / \ / / pn\ \ * / \ / / \ \ * / p1 \ p2 // 2 \\ * -------------- --------------- * p1 t0 p2 */ p.x=(Mx(tg->p0,M)+Mx(tg->p1,M)+Mx(tg->p2,M))*ONE_OVER_THREE ; p.y=(My(tg->p0,M)+My(tg->p1,M)+My(tg->p2,M))*ONE_OVER_THREE; if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) { fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } p.z[0]=(Mz(tg->p0,M)+Mz(tg->p1,M)+Mz(tg->p2,M))*ONE_OVER_THREE; p.val=VAR; point_number=BIT_insert(&M,p); if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles; tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1; tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2; tg1->boundary=tg2->boundary=tg3->boundary=tg->boundary; /* arrange neighbours */ tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2; tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0; tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1; /* new triangles are alive */ tg1->status=tg2->status=tg3->status=LIVE; /* old triangle is dead */ tg->status=DEAD; #ifdef DEBUG printf("subdivide_this: split to %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d)\n",tg1->n,tg1->p0,tg1->p1,tg1->p2,tg2->n,tg2->p0,tg2->p1,tg2->p2,tg3->n,tg3->p0,tg3->p1,tg3->p2); #endif /* update links from neighbours too */ tq.n=tg->t1;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg3->n); #endif update_neighbours_around(tempt,tg->n,tg3->n); tq.n=tg->t2;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg1->n); #endif update_neighbours_around(tempt,tg->n,tg1->n); tq.n=tg->t0;tempt=RBT_find(&tq, &rt); #ifdef DEBUG printf("subdivide_this: updating neighbour %d(%d,%d,%d) status %d as %d->%d\n",tempt->n,tempt->p0,tempt->p1,tempt->p2,tempt->status,tg->n,tg2->n); #endif update_neighbours_around(tempt,tg->n,tg2->n); parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2); tg2->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3); tg3->dag_p=current_dag_node; RBT_insert((void*)tg1, &rt); RBT_insert((void*)tg2, &rt); RBT_insert((void*)tg3, &rt); ntriangles+=3; /* since all edges are on boundary, we cannot refine them */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); delaunay_refine_constrained(tg3->n,point_number,&rt,&dt); return(1); } else if ( how_to_divide==4 ) { #ifdef DEBUG printf("subdivide_this: fixed_points=%d boundaries=%d badness="MDF" are="MDF"\n",fixed_points,boundary_edges,triangle_badness,area(tg)); #endif /* split the largest edge by drawing a perpendicular */ /* general form */ /* p0 p0 * /\ /\ * / \ t3 / |\ * t0 / \ / | \ * / t \ / 1| 2\ * p1 /____p___\ p3 p1/____|___\ p3 * \ / \ | 4 / * \ nei / \ 3| / * \ / \ | / * t1 \ / t2 \ |/ * \/ \/ * p2 p2 */ /* p1-p3 is the largest edge of triangle t */ /* if p1-p3 is a boundary, neighbour is hidden and ignored */ if ( (len2 >= len1) && (len2>=len0) ) { /* edge p0-p1 */ p0=tg->p2; p1=tg->p0; p3=tg->p1; tq.n=tg->t2; neigh=RBT_find(&tq,&rt); tq.n=tg->t1; t0=RBT_find(&tq,&rt); tq.n=tg->t0; t3=RBT_find(&tq,&rt); } else if ((len0 >= len1) && (len0>=len2) ) { /* edge p1-p2 */ p0=tg->p0; p1=tg->p1; p3=tg->p2; tq.n=tg->t0; neigh=RBT_find(&tq,&rt); tq.n=tg->t2; t0=RBT_find(&tq,&rt); tq.n=tg->t1; t3=RBT_find(&tq,&rt); } else { /* edge p2-p0 */ p0=tg->p1; p1=tg->p2; p3=tg->p0; tq.n=tg->t1; neigh=RBT_find(&tq,&rt); tq.n=tg->t0; t0=RBT_find(&tq,&rt); tq.n=tg->t2; t3=RBT_find(&tq,&rt); } /* determine foot of normal to longest edge */ /* note the denominator will be zero only if the triangle is * a right angled one. but then the longest edge will be the diagonal * hence in practice it will never be zero unless points coincide */ /* len0=-((Mx(p1)-Mx(p3))*(Mx(p3)-Mx(p0))+(My(p1)-My(p3))*(My(p3)-My(p0))) /((Mx(p1)-Mx(p0))*(Mx(p1)-Mx(p3))+(My(p1)-My(p0))*(My(p1)-My(p3))); p.x=(len0*Mx(p1)+Mx(p3))/(1+len0); p.y=(len0*My(p1)+My(p3))/(1+len0); */ p.x=(Mx(p1,M)+Mx(p3,M))*0.5; p.y=(My(p1,M)+My(p3,M))*0.5; if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) { fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__); exit(1); } /* p.z[0]=(len0*Mz(p1)+Mz(p3))/(1+len0); */ p.z[0]=(Mz(p1,M)+Mz(p3,M))*0.5; point_number=is_on_an_edge(p1,p3); /* point_number= boundary_edge+1, if (p1,p3) lies on an edge */ p.val=(point_number && edge_array[point_number-1].type==DR?FX:VAR); #ifdef DEBUG if (p.val==FX) { printf("subdivide_this: point fixed\n"); } else { printf("subdivide_this: point variable\n"); } #endif point_number=BIT_insert(&M,p); /* check for neighbour */ if (neigh && neigh->status==LIVE) { if (neigh->t0==tg->n) { p2=neigh->p0; tq.n=neigh->t1; t1=RBT_find(&tq,&rt); tq.n=neigh->t2; t2=RBT_find(&tq,&rt); } else if (neigh->t1 == tg->n) { p2=neigh->p1; tq.n=neigh->t2; t1=RBT_find(&tq,&rt); tq.n=neigh->t0; t2=RBT_find(&tq,&rt); } else { p2=neigh->p2; tq.n=neigh->t0; t1=RBT_find(&tq,&rt); tq.n=neigh->t1; t2=RBT_find(&tq,&rt); } #ifdef DEBUG printf("subdivide_this: neighbour %d with neighbours (%d(%d),%d(%d),%d(%d))\n",neigh->n,neigh->t0,neigh->p0,neigh->t1,neigh->p1,neigh->t2,neigh->p2); printf("sundivide_this: generalized to p0=%d,p1=%d,p2=%d,p3=%d\n",p0,p1,p2,p3); #endif } /*split to 4 or 2 */ if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ( neigh && neigh->status==LIVE) { if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) { fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__); exit(1); } } /* update triangle data */ tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles; tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1; tg1->boundary=tg2->boundary=tg->boundary; /* alive, dead triangles */ tg1->status=tg2->status=LIVE; tg->status=DEAD; if (neigh && neigh->status==LIVE) { tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2; tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3; tg3->status=tg4->status=LIVE; tg3->boundary=tg4->boundary=neigh->boundary; /* at the last moment this */ neigh->status=DEAD; tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1); tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n; tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n; tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n; update_neighbours_around(t0,tg->n,tg1->n); update_neighbours_around(t3,tg->n,tg2->n); update_neighbours_around(t1,neigh->n,tg3->n); update_neighbours_around(t2,neigh->n,tg4->n); /* insert triangles into RBT */ RBT_insert((void*)tg1,&rt); RBT_insert((void*)tg2,&rt); RBT_insert((void*)tg3,&rt); RBT_insert((void*)tg4,&rt); ntriangles+=4; /* update DAG */ parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2); tg2->dag_p=current_dag_node; parent_dag_node=neigh->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg3); tg3->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg4); tg4->dag_p=current_dag_node; #ifdef DEBUG printf("split to 4 triangles %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d), %d(%d,%d,%d)\n",tg1->n,tg1->p0,tg1->p1,tg1->p2,tg2->n,tg2->p0,tg2->p1,tg2->p2,tg3->n,tg3->p0,tg3->p1,tg3->p2,tg4->n,tg4->p0,tg4->p1,tg4->p2); #endif /* refine the triangulation */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); delaunay_refine_constrained(tg3->n,point_number,&rt,&dt); delaunay_refine_constrained(tg4->n,point_number,&rt,&dt); } else { tg1->t0=(neigh?neigh->n:-1);tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1); tg2->t0=(neigh?neigh->n:-1);tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n; update_neighbours_around(t0,tg->n,tg1->n); update_neighbours_around(t3,tg->n,tg2->n); /* insert triangles into RBT */ RBT_insert((void*)tg1,&rt); RBT_insert((void*)tg2,&rt); ntriangles+=2; /* update DAG */ parent_dag_node=tg->dag_p; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1); tg1->dag_p=current_dag_node; current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2); tg2->dag_p=current_dag_node; /* refine the triangulation */ delaunay_refine_constrained(tg1->n,point_number,&rt,&dt); delaunay_refine_constrained(tg2->n,point_number,&rt,&dt); } return(1); } return(0); } /* similar to subdivide_all_triangles * but we take into account the gradient * at each point when dividing a triangle */ int refine_mesh_with_iterations(MY_DOUBLE priority) { triangle *rec; int flag=0; /* gradient parameters */ MY_DOUBLE grad_mean,max_grad,triangle_gradient; MY_INT count; /* first do a traversal to see the max-min gradients are */ grad_mean=max_grad=0; count=0; DAG_reverse_traverse_list_reset(&dt); rec=DAG_reverse_traverse_prune_list(&dt); while(rec){ if (rec &&(rec->status==LIVE) ) { triangle_gradient=ABS(Mz(rec->p0,M)-Mz(rec->p1,M)) + ABS(Mz(rec->p1,M)-Mz(rec->p2,M)) + ABS(Mz(rec->p2,M)-Mz(rec->p0,M)); if ( max_grad < triangle_gradient ) max_grad=triangle_gradient; grad_mean=((MY_DOUBLE)count*grad_mean+triangle_gradient)/((MY_DOUBLE)count+1.0); count++; #ifdef DEBUG printf("refine_mesh_iter: %d gradient"MDF"\n",rec->n,triangle_gradient); #endif } rec=DAG_reverse_traverse_prune_list(&dt); } /* refine all triangles with gradient >1/2(mean+max) */ g_gradient_limit=0.5*(max_grad+grad_mean); #ifdef DEBUG printf("refine_mesh_iter: gradient max,mean "MDF"and "MDF"\n",max_grad,grad_mean); printf("refine_mesh_iter: set gradient limit"MDF"\n",g_gradient_limit); #endif /* now do the proper refinement */ DAG_reverse_traverse_list_reset(&dt); rec=DAG_reverse_traverse_prune_list(&dt); while(rec){ if (rec &&(rec->status==LIVE) && refine_this_triangle(rec,0,priority)) { /* if a triangle was divided, flag=1 */ #ifdef DEBUG printf("refining %d\n",rec->n); #endif flag |=1; } rec=DAG_reverse_traverse_prune_list(&dt); } return(flag); }