1 /* pdnmesh - a 2D finite element solver
2 Copyright (C) 2001-2005 Sarod Yatawatta <sarod@users.sf.net>
3 This program is free software; you can redistribute it and/or modify
4 it under the terms of the GNU General Public License as published by
5 the Free Software Foundation; either version 2 of the License, or
6 (at your option) any later version.
7
8 This program is distributed in the hope that it will be useful,
9 but WITHOUT ANY WARRANTY; without even the implied warranty of
10 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11 GNU General Public License for more details.
12
13 You should have received a copy of the GNU General Public License
14 along with this program; if not, write to the Free Software
15 Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
16 $Id: subdivide.c,v 1.25 2005/04/18 08:30:00 sarod Exp $
17 */
18
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23
24 /* subdivision and refinement of the mesh */
25 #include <stdio.h>
26 #include <stdlib.h>
27 #include "types.h"
28
29
30 /* parameters for modifying the mesh */
31 MY_DOUBLE g_badness_limit;
32 MY_DOUBLE g_area_floor;
33 MY_DOUBLE g_area_ceil;
34 MY_DOUBLE g_gradient_limit;
35
36 /* determines line defined by (p1,p2) is
37 * on an edge */
38 /* returns 0 if no edge, (edge no. +1) if is on an edge */
39 static int
is_on_an_edge(MY_INT p1,MY_INT p2)40 is_on_an_edge(MY_INT p1,MY_INT p2)
41 {
42 MY_INT i;
43 MY_DOUBLE x1,x2,y1,y2;
44 for (i=0; i<nedges; i++){
45 x1=Mx(edge_array[i].p1,M);
46 x2=Mx(edge_array[i].p2,M);
47 y1=My(edge_array[i].p1,M);
48 y2=My(edge_array[i].p2,M);
49 #ifdef DEBUG
50 printf("is_on_an_edge: %d with ("MDF","MDF")-("MDF","MDF") on ? %d("MDF","MDF"),%d("MDF","MDF")\n",i,x1,y1,x2,y2,p1,Mx(p1,M),My(p1,M),p2,Mx(p2,M),My(p2,M));
51 #endif
52 /* first consider coincident points case */
53 if ( (IS_ZERO(Mx(p1,M)-x1) && IS_ZERO(My(p1,M)-y1)
54 && IS_ZERO(Mx(p2,M)-x2) && IS_ZERO(My(p2,M)-y2) )
55 || (IS_ZERO(Mx(p1,M)-x2) && IS_ZERO(My(p1,M)-y2)
56 && IS_ZERO(Mx(p2,M)-x1) && IS_ZERO(My(p2,M)-y1) ) ) {
57 #ifdef DEBUG
58 printf(".. Yes\n");
59 #endif
60 return(i+1);
61 }
62
63 #ifdef DEBUG
64 printf(MDF", "MDF"::",((Mx(p1,M)-x1)*(y1-y2)-(My(p1,M)-y1)*(x1-x2)),((Mx(p2,M)-x1)*(y1-y2)-(My(p2,M)-y1)*(x1-x2)));
65 printf("%d, %d\n",IS_ZERO((Mx(p1,M)-x1)*(y1-y2)-(My(p1,M)-y1)*(x1-x2)),IS_ZERO((Mx(p2,M)-x1)*(y1-y2)-(My(p2,M)-y1)*(x1-x2)));
66 #endif
67 if ( IS_ZERO((Mx(p1,M)-x1)*(y1-y2)-(My(p1,M)-y1)*(x1-x2))
68 && IS_ZERO((Mx(p2,M)-x1)*(y1-y2)-(My(p2,M)-y1)*(x1-x2)) ) {
69 /* if we are here, the line is parallel to an edge */
70 /* now both points must be within the edge */
71 #ifdef DEBUG
72 printf(MDF" "MDF" "MDF" "MDF"\n",(Mx(p1,M)-x2)*(x1-Mx(p1,M)),(My(p1,M)-y2)*(y1-My(p1,M)),(Mx(p2,M)-x2)*(x1-Mx(p2,M)),(My(p2,M)-y2)*(y1-My(p2,M)));
73 #endif
74 if ( ((Mx(p1,M)-x2)*(x1-Mx(p1,M))>= 0)
75 && ((My(p1,M)-y2)*(y1-My(p1,M))>= 0)
76 && ((Mx(p2,M)-x2)*(x1-Mx(p2,M))>= 0)
77 && ((My(p2,M)-y2)*(y1-My(p2,M))>= 0) ) {
78 #ifdef DEBUG
79 printf(".. Yes\n");
80 #endif
81 return(i+1);
82 }
83 }
84 }
85 #ifdef DEBUG
86 printf(".. NO\n");
87 #endif
88
89 return(0);
90 }
91
92
93 /* a function to calculate the area of a triangle */
94 static MY_DOUBLE
area(triangle * t)95 area( triangle *t )
96 {
97 MY_DOUBLE temp1;
98 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));
99 return(0.5*ABS(temp1));
100 }
101
102 /* update neighbour link in a triangle */
103 /* neighbour old will be replaced by neighbour new */
104 static void
update_neighbours_around(triangle * neighbour,int old,int new)105 update_neighbours_around(triangle *neighbour,int old,int new)
106 {
107 if ( neighbour ) {
108 if ( neighbour->t0 == old ) {
109 neighbour->t0 = new;
110 } else if ( neighbour->t1 == old ) {
111 neighbour->t1 = new;
112 } else if ( neighbour->t2 == old ) {
113 neighbour->t2 = new;
114 }
115 }
116 }
117
118 /* a function to determine whether point p4 is inside */
119 /* the circumcircle of points p1,p2,p3 */
120 static int
inside_circle(int p1,int p2,int p3,int p4)121 inside_circle( int p1, int p2, int p3, int p4)
122 {
123 double temp1,temp2,temp3,temp4,a2,b2;
124 #ifdef DEBUG
125 printf("checking (%d,%d,%d) includes %d\n",p1,p2,p3,p4);
126 #endif
127 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));
128 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));
129 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));
130 /* we know this is not zero */
131 temp3 = 1/temp3;
132 a2 = (temp1*(My(p1,M)-My(p3,M))-temp2*(My(p1,M)-My(p2,M)))*temp3;
133 b2 = (temp2*(Mx(p1,M)-Mx(p2,M))-temp1*(Mx(p1,M)-Mx(p3,M)))*temp3;
134 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);
135 if ( a2*(Mx(p1,M)-Mx(p4,M))+b2*(My(p1,M)-My(p4,M)) < temp4 ) {
136 #ifdef DEBUG
137 printf("points are inside circle\n");
138 #endif
139 return(1);
140 }
141 return(0);
142 }
143
144 /* tgn - triangle number to refine */
145 /* ptn - point number of vertex of tgn facing the edge to be refined */
146 /* assertion - at each step, the remainder of triangles form a valid
147 * delaunay triangulation */
148 int
delaunay_refine_constrained(int tgn,int ptn,RBT_tree * rbt,DAG_tree * dag)149 delaunay_refine_constrained(int tgn, int ptn, RBT_tree * rbt,DAG_tree *dag)
150 {
151 triangle tq,*t,*t_n;
152 int neighbour,p0,p1,p2,p4;
153 int t1,t2,t3,t4; /* neighbours */
154 triangle *tp1,*tp2,*tp3,*tp4,*tnew1,*tnew2;
155
156 DAG_node *current_dag_node,*parent_dag_node;
157 /* fix everything into this generalized form
158 * p2/\
159 * /| \
160 * / | \
161 * t1/ | \ t4
162 * / | \
163 * p0 / | \ p4
164 * \ t | tn /
165 * \ | /
166 * t2 \ | / t3
167 * \ | /
168 * p1
169 */
170
171 /* get the triangle */
172 tq.n=tgn;
173 t=RBT_find(&tq, rbt);
174 if ( t && ( t->n == tgn )) { /* triangle exists */
175 /* first find neighbour facing the edge to be refined */
176 if ( ptn == t->p0 ) {
177 neighbour=t->t0;
178 p0=t->p0;p1=t->p1;p2=t->p2;
179 t1=t->t1; t2=t->t2;
180 } else if ( ptn == t->p1 ) {
181 neighbour=t->t1;
182 p0=t->p1;p1=t->p2;p2=t->p0;
183 t1=t->t2; t2=t->t0;
184 } else {
185 neighbour=t->t2;
186 p0=t->p2;p1=t->p0;p2=t->p1;
187 t1=t->t0; t2=t->t1;
188 }
189 /* if the edge to be flipped is on a boundary,
190 * we cannot do anything, we quit */
191 if ( is_on_an_edge(p1,p2) ) {
192 #ifdef DEBUG
193 printf("delaunay_refine_cons: edge on boundary :cannot refine\n");
194 #endif
195 return(1);
196 }
197 /* now find the neighbouring triangle */
198 /* it may not exist */
199 tq.n=neighbour;
200 t_n=RBT_find(&tq, rbt);
201 if ( (!t_n) || (t_n->status!=LIVE) ) {/* no neighbour */
202 return(0);
203 }
204 /* find the shared edge with neigbour, or rather, the opposing vertex */
205 if ( tgn == t_n->t0 ) {
206 p4=t_n->p0;
207 t3=t_n->t1; t4=t_n->t2;
208 } else if ( tgn == t_n->t1 ) {
209 p4=t_n->p1;
210 t3=t_n->t2; t4=t_n->t0;
211 } else {
212 p4=t_n->p2;
213 t3=t_n->t0; t4=t_n->t1;
214 }
215
216 /* if we form a triangle with all three vertices fixed, we cannot refine */
217 if ((Mval(p0,M)==FX)&&(Mval(p1,M)==FX)&&(Mval(p4,M)==FX)) {
218 #ifdef DEBUG
219 printf("delaunay_refine_cons: fixed points: cannot refine\n");
220 #endif
221 return(1);
222 }
223
224 /* now we have fitted everything to the general framework */
225 if ( inside_circle(p0,p2,p1,p4) ) {
226 /* we flip the edge (p1,p2) */
227 /* this could only be done if p4 can be seen
228 * by p0, i.e. (p1,p2) is the largest edge of neighbour
229 * but since we assert a valid initial triangulation,
230 * this will hold.
231 */
232
233 /* get neighbouring triangles first */
234 tq.n=t1;tp1=RBT_find(&tq, rbt);
235 tq.n=t2;tp2=RBT_find(&tq, rbt);
236 tq.n=t3;tp3=RBT_find(&tq, rbt);
237 tq.n=t4;tp4=RBT_find(&tq, rbt);
238 /* create 2 new triangles */
239 if ((tnew1=(triangle *)malloc(sizeof(triangle)))==0) {
240 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
241 exit(1);
242 }
243 if ((tnew2=(triangle *)malloc(sizeof(triangle)))==0) {
244 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
245 exit(1);
246 }
247 /* we transform the setup into
248 * p2/\
249 * / \
250 * / \
251 * t1/ \ t4
252 * / new1 \
253 * p0 / \ p4
254 * \- - ------ /
255 * \ new2 /
256 * t2 \ / t3
257 * \ /
258 * p1
259 */
260 tnew1->p0=p0;tnew1->p1=p4;tnew1->p2=p2;
261 tnew1->n=ntriangles;
262 tnew2->p0=p0;tnew2->p1=p1;tnew2->p2=p4;
263 tnew2->n=ntriangles+1;
264 tnew1->status=tnew2->status=LIVE;
265 tnew1->boundary=tnew2->boundary=t->boundary;
266 t->status=t_n->status=DEAD;
267 /* neighbours */
268 tnew1->t0=t4;tnew1->t1=t1;tnew1->t2=tnew2->n;
269 tnew2->t0=t3;tnew2->t1=tnew1->n;tnew2->t2=t2;
270 /* undate the neighbours of the neighbours too */
271 update_neighbours_around(tp1,tgn,tnew1->n);
272 update_neighbours_around(tp2,tgn,tnew2->n);
273 update_neighbours_around(tp3,neighbour,tnew2->n);
274 update_neighbours_around(tp4,neighbour,tnew1->n);
275 /* insert the 2 new triangles into RB tree */
276 RBT_insert((void*)tnew1, rbt);
277
278 parent_dag_node=t->dag_p;
279 current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew1);
280 tnew1->dag_p=current_dag_node;
281 /* link it to the other parent in dag */
282 DAG_link(dag,t_n->dag_p,current_dag_node);
283
284 /* second triangle */
285 RBT_insert((void*)tnew2, rbt);
286
287 parent_dag_node=t->dag_p;
288 current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew2);
289 tnew2->dag_p=current_dag_node;
290 /* link it to the other parent in dag */
291 DAG_link(dag,t_n->dag_p,current_dag_node);
292
293 /* increment triangle count */
294 ntriangles +=2;
295 /* printf("flipped %d neighbour %d to %d and %d\n",tgn,neighbour,tnew1->n,tnew2->n); */
296
297 /* call this function recursively */
298 delaunay_refine_constrained(tnew1->n,p0,rbt,dag);
299 delaunay_refine_constrained(tnew2->n,p0,rbt,dag);
300 }
301
302
303 } else {
304 /* error */
305 printf("delaunay_refine_cons: cannot find triangle. something is wrong\n");
306 return(-1);
307 }
308
309
310 return(0);
311 }
312
313 /* decide to split triangle or not */
314 static int
do_we_split_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness)315 do_we_split_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness)
316 {
317 #ifdef DEBUG
318 printf("do_we_split: badness "MDF">? "MDF", "MDF" > area > "MDF"\n",badness, g_badness_limit,g_area_floor,g_area_ceil);
319 #endif
320 if (( (badness>g_badness_limit) && (tarea>g_area_floor))
321 || (tarea >g_area_ceil) ) {
322 return(1);
323 }
324 return(0);
325 }
326
327 /* decide to split triangle or not */
328 /* also consider gradient */
329 static int
do_we_refine_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness,MY_DOUBLE gradient)330 do_we_refine_this_triangle(MY_DOUBLE tarea,MY_DOUBLE badness,MY_DOUBLE gradient)
331 {
332 #ifdef DEBUG
333 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);
334 #endif
335 if ( gradient > g_gradient_limit
336 || ((badness>g_badness_limit) && (tarea>g_area_floor))
337 || (tarea >g_area_ceil) ) {
338 return(1);
339 }
340 return(0);
341 }
342 /* given triangle t and its neighbour n
343 * determine whether we can split to 4. we can when
344 * they share a boundary edge or the edge
345 * common is the largest of both triangles
346 */
347 static int
can_we_split_to_4(triangle * t,triangle * n)348 can_we_split_to_4(triangle *t, triangle *n) {
349 if(t->boundary != n->boundary ) {
350 return(1);
351 }
352 /* find the largest edges */
353 /* we know the largest edge of t is one common with n */
354 /* find the largest edge of n */
355 if ( n->t0==t->n ) {
356 return((LENGTH(n->p1,n->p2,M)>=LENGTH(n->p0,n->p2,M))
357 && (LENGTH(n->p1,n->p2,M)>=LENGTH(n->p0,n->p1,M)));
358 }
359 if ( n->t1==t->n ) {
360 return((LENGTH(n->p2,n->p0,M)>=LENGTH(n->p1,n->p2,M))
361 && (LENGTH(n->p2,n->p0,M)>=LENGTH(n->p1,n->p0,M)));
362 }
363 /* else n->t2==t->n */
364 return((LENGTH(n->p0,n->p1,M)>=LENGTH(n->p1,n->p2,M))
365 && (LENGTH(n->p0,n->p1,M)>=LENGTH(n->p0,n->p2,M)));
366 }
367
368
369 /* subdivision of a triangle */
370 int
subdivide_this_triangle(triangle * tg)371 subdivide_this_triangle(triangle *tg)
372 {
373 /* our policy is simple and as follows:
374 * 1) all edges on boundary - insert point at centroid - split to 3
375 * 2) in all other cases, draw a perpendicular to one of the edges
376 * split to 4 or 2
377 * the perpendicular is drawn to the longest edge, if other two edges
378 * do not belong to the boundary, else, the free edge is chosen to draw it
379 */
380 int fixed_points, boundary_edges;
381 point p;
382 DAG_node *current_dag_node,*parent_dag_node;
383 triangle *tg1,*tg2,*tg3,*tg4,*tempt,tq;
384 triangle *t0,*t1,*t2,*t3,*neigh;
385 MY_INT point_number,p0,p1,p2,p3;
386 MY_DOUBLE len0,len1,len2;
387 MY_DOUBLE min_edge,max_edge,triangle_badness;
388
389 /* how to divide ? */
390 int how_to_divide;
391 /* 3: divide into 3, inserting point at centroid
392 * 4: divide into 4 (or 2) inserting a normal to longest edge
393 * 0: no division */
394
395 /* ignore some warnings about uninitialized variables */
396 t1=t2=tg3=tg4=0;
397 p2=0;
398
399 fixed_points=0;
400 if (Mval(tg->p0,M)==FX) {
401 fixed_points++;
402 }
403 if (Mval(tg->p1,M)==FX) {
404 fixed_points++;
405 }
406 if (Mval(tg->p2,M)==FX) {
407 fixed_points++;
408 }
409 boundary_edges=0;
410 if (is_on_an_edge(tg->p0,tg->p1)) {
411 boundary_edges++;
412 }
413 if (is_on_an_edge(tg->p1,tg->p2)) {
414 boundary_edges++;
415 }
416 if (is_on_an_edge(tg->p2,tg->p0)) {
417 boundary_edges++;
418 }
419 /* lengths of edges */
420 len0=LENGTH(tg->p1,tg->p2,M);
421 len1=LENGTH(tg->p2,tg->p0,M);
422 len2=LENGTH(tg->p1,tg->p0,M);
423 min_edge=(((len0<=len1) && (len0<=len1))? len0 : \
424 (((len1<=len0) && (len1<=len2))? len1 : len2));
425 max_edge=(((len0>=len1) && (len0>=len1))? len0 : \
426 (((len1>=len0) && (len1>=len2))? len1 : len2));
427
428 triangle_badness=max_edge/min_edge;
429 #ifdef DEBUG
430 printf("\n\nsubdivide_this: considering triangle %d points(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2);
431 printf("subdivide_this: min edge "MDF", max edge "MDF", badness "MDF" area "MDF"\n",min_edge,max_edge,triangle_badness, area(tg));
432 printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges);
433 #endif
434 how_to_divide=0;
435 if ( fixed_points == 3 ) {
436 how_to_divide=3;
437 } else if ( do_we_split_this_triangle(area(tg),triangle_badness) ) {
438 /* find neighbour */
439 if ( (len2 >= len1) && (len2>=len0) ) {
440 tq.n=tg->t2;
441 neigh=RBT_find(&tq,&rt);
442
443 } else if ((len0 >= len1) && (len0>=len2) ) {
444 tq.n=tg->t0;
445 neigh=RBT_find(&tq,&rt);
446
447 } else {
448 tq.n=tg->t1;
449 neigh=RBT_find(&tq,&rt);
450 }
451 if ( (!neigh) || (neigh->status!=LIVE ) ) {
452 /* no neighbour exists, divide into 4 or 2 */
453 how_to_divide=4;
454 } else if ( neigh &&can_we_split_to_4(tg, neigh) ) {
455 /* if the largest edge is a boundary edge
456 * we know it cannot be flipped by delaunay refinement.
457 * so we split to 4
458 * also if the larges edge of neighbour is the largest of triangle,
459 * we split to 4
460 */
461 how_to_divide=4;
462 } else {
463 how_to_divide=3;
464 }
465 /* if the neighbours largest edge is also our largest edge,
466 * we can split to 4 */
467 }
468
469 if ( how_to_divide==3 ) {
470 #ifdef DEBUG
471 printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges);
472 #endif
473 /* split to 3 by inserting point at centroid */
474 /* p0 /\ p0 /|\
475 * / \ / | \
476 * / \ to /1 | \ t1
477 * / \ t2 / / \3 \
478 * / \ / / pn\ \
479 * / \ / / \ \
480 * / p1 \ p2 // 2 \\
481 * -------------- ---------------
482 * p1 t0 p2
483 */
484 p.x=(Mx(tg->p0,M)+Mx(tg->p1,M)+Mx(tg->p2,M))*ONE_OVER_THREE;
485 p.y=(My(tg->p0,M)+My(tg->p1,M)+My(tg->p2,M))*ONE_OVER_THREE;
486 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
487 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
488 exit(1);
489 }
490 p.z[0]=(Mz(tg->p0,M)+Mz(tg->p1,M)+Mz(tg->p2,M))*ONE_OVER_THREE;
491 p.val=VAR;
492
493 point_number=BIT_insert(&M,p);
494 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
495 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
496 exit(1);
497 }
498 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
499 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
500 exit(1);
501 }
502 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
503 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
504 exit(1);
505 }
506 tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
507 tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
508 tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
509 tg1->boundary=tg2->boundary=tg3->boundary=tg->boundary;
510 /* arrange neighbours */
511 tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
512 tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
513 tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
514 /* new triangles are alive */
515 tg1->status=tg2->status=tg3->status=LIVE;
516 /* old triangle is dead */
517 tg->status=DEAD;
518 #ifdef DEBUG
519 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);
520 #endif
521
522 /* update links from neighbours too */
523 tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
524 #ifdef DEBUG
525 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);
526 #endif
527 update_neighbours_around(tempt,tg->n,tg3->n);
528 tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
529 #ifdef DEBUG
530 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);
531 #endif
532 update_neighbours_around(tempt,tg->n,tg1->n);
533 tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
534 #ifdef DEBUG
535 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);
536 #endif
537 update_neighbours_around(tempt,tg->n,tg2->n);
538
539
540 parent_dag_node=tg->dag_p;
541 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
542 tg1->dag_p=current_dag_node;
543 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
544 tg2->dag_p=current_dag_node;
545 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
546 tg3->dag_p=current_dag_node;
547
548 RBT_insert((void*)tg1, &rt);
549 RBT_insert((void*)tg2, &rt);
550 RBT_insert((void*)tg3, &rt);
551 ntriangles+=3;
552
553 /* since all edges are on boundary, we cannot refine them */
554 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
555 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
556 delaunay_refine_constrained(tg3->n,point_number,&rt,&dt);
557
558 return(1);
559
560 } else if ( how_to_divide==4 ) {
561 #ifdef DEBUG
562 printf("subdivide_this: fixed_points=%d boundaries=%d badness="MDF" are="MDF"\n",fixed_points,boundary_edges,triangle_badness,area(tg));
563 #endif
564
565 /* split the largest edge by drawing a perpendicular */
566 /* general form */
567 /* p0 p0
568 * /\ /\
569 * / \ t3 / |\
570 * t0 / \ / | \
571 * / t \ / 1| 2\
572 * p1 /____p___\ p3 p1/____|___\ p3
573 * \ / \ | 4 /
574 * \ nei / \ 3| /
575 * \ / \ | /
576 * t1 \ / t2 \ |/
577 * \/ \/
578 * p2 p2
579 */
580 /* p1-p3 is the largest edge of triangle t */
581 /* if p1-p3 is a boundary, neighbour is hidden and ignored */
582 if ( (len2 >= len1) && (len2>=len0) ) {
583 /* edge p0-p1 */
584 p0=tg->p2;
585 p1=tg->p0;
586 p3=tg->p1;
587 tq.n=tg->t2;
588 neigh=RBT_find(&tq,&rt);
589
590 tq.n=tg->t1;
591 t0=RBT_find(&tq,&rt);
592 tq.n=tg->t0;
593 t3=RBT_find(&tq,&rt);
594 } else if ((len0 >= len1) && (len0>=len2) ) {
595 /* edge p1-p2 */
596 p0=tg->p0;
597 p1=tg->p1;
598 p3=tg->p2;
599 tq.n=tg->t0;
600 neigh=RBT_find(&tq,&rt);
601
602 tq.n=tg->t2;
603 t0=RBT_find(&tq,&rt);
604 tq.n=tg->t1;
605 t3=RBT_find(&tq,&rt);
606 } else {
607 /* edge p2-p0 */
608 p0=tg->p1;
609 p1=tg->p2;
610 p3=tg->p0;
611 tq.n=tg->t1;
612 neigh=RBT_find(&tq,&rt);
613
614 tq.n=tg->t0;
615 t0=RBT_find(&tq,&rt);
616 tq.n=tg->t2;
617 t3=RBT_find(&tq,&rt);
618 }
619 /* determine foot of normal to longest edge */
620 /* note the denominator will be zero only if the triangle is
621 * a right angled one. but then the longest edge will be the diagonal
622 * hence in practice it will never be zero unless points coincide */
623 /* len0=-((Mx(p1)-Mx(p3))*(Mx(p3)-Mx(p0))+(My(p1)-My(p3))*(My(p3)-My(p0)))
624 /((Mx(p1)-Mx(p0))*(Mx(p1)-Mx(p3))+(My(p1)-My(p0))*(My(p1)-My(p3)));
625 p.x=(len0*Mx(p1)+Mx(p3))/(1+len0);
626 p.y=(len0*My(p1)+My(p3))/(1+len0); */
627 p.x=(Mx(p1,M)+Mx(p3,M))*0.5;
628 p.y=(My(p1,M)+My(p3,M))*0.5;
629 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
630 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
631 exit(1);
632 }
633 /* p.z[0]=(len0*Mz(p1)+Mz(p3))/(1+len0); */
634 p.z[0]=(Mz(p1,M)+Mz(p3,M))*0.5;
635 point_number=is_on_an_edge(p1,p3);
636 /* point_number= boundary_edge+1, if (p1,p3) lies on an edge */
637 p.val=(point_number && edge_array[point_number-1].type==DR?FX:VAR);
638 #ifdef DEBUG
639 if (p.val==FX) {
640 printf("subdivide_this: point fixed\n");
641 } else {
642 printf("subdivide_this: point variable\n");
643 }
644 #endif
645 point_number=BIT_insert(&M,p);
646 /* check for neighbour */
647 if (neigh && neigh->status==LIVE) {
648 if (neigh->t0==tg->n) {
649 p2=neigh->p0;
650 tq.n=neigh->t1;
651 t1=RBT_find(&tq,&rt);
652 tq.n=neigh->t2;
653 t2=RBT_find(&tq,&rt);
654 } else if (neigh->t1 == tg->n) {
655 p2=neigh->p1;
656 tq.n=neigh->t2;
657 t1=RBT_find(&tq,&rt);
658 tq.n=neigh->t0;
659 t2=RBT_find(&tq,&rt);
660 } else {
661 p2=neigh->p2;
662 tq.n=neigh->t0;
663 t1=RBT_find(&tq,&rt);
664 tq.n=neigh->t1;
665 t2=RBT_find(&tq,&rt);
666 }
667 #ifdef DEBUG
668 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);
669 printf("sundivide_this: generalized to p0=%d,p1=%d,p2=%d,p3=%d\n",p0,p1,p2,p3);
670 #endif
671 }
672 /*split to 4 or 2 */
673 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
674 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
675 exit(1);
676 }
677 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
678 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
679 exit(1);
680 }
681 if ( neigh && neigh->status==LIVE) {
682 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
683 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
684 exit(1);
685 }
686 if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
687 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
688 exit(1);
689 }
690 }
691
692 /* update triangle data */
693 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
694 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
695 tg1->boundary=tg2->boundary=tg->boundary;
696 /* alive, dead triangles */
697 tg1->status=tg2->status=LIVE;
698 tg->status=DEAD;
699 if (neigh && neigh->status==LIVE) {
700 tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
701 tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;
702 tg3->status=tg4->status=LIVE;
703 tg3->boundary=tg4->boundary=neigh->boundary;
704 /* at the last moment this */
705 neigh->status=DEAD;
706 tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
707 tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
708 tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
709 tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
710
711 update_neighbours_around(t0,tg->n,tg1->n);
712 update_neighbours_around(t3,tg->n,tg2->n);
713 update_neighbours_around(t1,neigh->n,tg3->n);
714 update_neighbours_around(t2,neigh->n,tg4->n);
715 /* insert triangles into RBT */
716 RBT_insert((void*)tg1,&rt);
717 RBT_insert((void*)tg2,&rt);
718 RBT_insert((void*)tg3,&rt);
719 RBT_insert((void*)tg4,&rt);
720 ntriangles+=4;
721 /* update DAG */
722 parent_dag_node=tg->dag_p;
723 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1);
724 tg1->dag_p=current_dag_node;
725 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2);
726 tg2->dag_p=current_dag_node;
727 parent_dag_node=neigh->dag_p;
728 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg3);
729 tg3->dag_p=current_dag_node;
730 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg4);
731 tg4->dag_p=current_dag_node;
732 #ifdef DEBUG
733 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);
734 #endif
735 /* refine the triangulation */
736 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
737 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
738 delaunay_refine_constrained(tg3->n,point_number,&rt,&dt);
739 delaunay_refine_constrained(tg4->n,point_number,&rt,&dt);
740
741
742 } else {
743 tg1->t0=(neigh?neigh->n:-1);tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
744 tg2->t0=(neigh?neigh->n:-1);tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
745 update_neighbours_around(t0,tg->n,tg1->n);
746 update_neighbours_around(t3,tg->n,tg2->n);
747 /* insert triangles into RBT */
748 RBT_insert((void*)tg1,&rt);
749 RBT_insert((void*)tg2,&rt);
750 ntriangles+=2;
751 /* update DAG */
752 parent_dag_node=tg->dag_p;
753 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1);
754 tg1->dag_p=current_dag_node;
755 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2);
756 tg2->dag_p=current_dag_node;
757
758 /* refine the triangulation */
759 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
760 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
761
762 }
763 return(1);
764
765 }
766
767 return(0);
768 }
769
770 /* subdivides all live triangles */
771 int
subdivide_all_triangles(void)772 subdivide_all_triangles(void)
773 {
774 triangle *rec;
775 DAG_traverse_list_reset(&dt);
776 rec=DAG_traverse_prune_list(&dt);
777 while(rec){
778 if (rec &&(rec->status==LIVE)) {
779 subdivide_this_triangle(rec);
780 }
781 rec=DAG_traverse_prune_list(&dt);
782 }
783 return(0);
784 }
785
786 /* refinement of a triangle */
787 /* similar to sundivide_this_triangle
788 * except we consired gradient of potentials */
789 /* if override==1, we split it no matter what */
790 int
refine_this_triangle(triangle * tg,int override,MY_DOUBLE priority)791 refine_this_triangle(triangle *tg, int override,MY_DOUBLE priority)
792 {
793 /* our policy is simple and as follows:
794 * 1) all edges on boundary - insert point at centroid - split to 3
795 * 2) in all other cases, draw a perpendicular to one of the edges
796 * split to 4 or 2
797 * the perpendicular is drawn to the longest edge, if other two edges
798 * do not belong to the boundary, else, the free edge is chosen to draw it
799 */
800 int fixed_points, boundary_edges;
801 point p;
802 DAG_node *current_dag_node,*parent_dag_node;
803 triangle *tg1,*tg2,*tg3,*tg4,*tempt,tq;
804 triangle *t0,*t1,*t2,*t3,*neigh;
805 MY_INT point_number,p0,p1,p2,p3;
806 MY_DOUBLE len0,len1,len2;
807 MY_DOUBLE min_edge,max_edge,triangle_badness,triangle_gradient;
808
809 /* how to divide ? */
810 int how_to_divide;
811 /* 3: divide into 3, inserting point at centroid
812 * 4: divide into 4 (or 2) inserting a normal to longest edge
813 * 0: no division */
814
815 /* ignore some warnings about uninitialized variables */
816 t1=t2=tg3=tg4=0;
817 p2=0;
818
819 fixed_points=0;
820 if (Mval(tg->p0,M)==FX) {
821 fixed_points++;
822 }
823 if (Mval(tg->p1,M)==FX) {
824 fixed_points++;
825 }
826 if (Mval(tg->p2,M)==FX) {
827 fixed_points++;
828 }
829 boundary_edges=0;
830 if (is_on_an_edge(tg->p0,tg->p1)) {
831 boundary_edges++;
832 }
833 if (is_on_an_edge(tg->p1,tg->p2)) {
834 boundary_edges++;
835 }
836 if (is_on_an_edge(tg->p2,tg->p0)) {
837 boundary_edges++;
838 }
839 /* lengths of edges */
840 len0=LENGTH(tg->p1,tg->p2,M);
841 len1=LENGTH(tg->p2,tg->p0,M);
842 len2=LENGTH(tg->p1,tg->p0,M);
843 min_edge=(((len0<=len1) && (len0<=len1))? len0 : \
844 (((len1<=len0) && (len1<=len2))? len1 : len2));
845 max_edge=(((len0>=len1) && (len0>=len1))? len0 : \
846 (((len1>=len0) && (len1>=len2))? len1 : len2));
847
848 triangle_badness=max_edge/min_edge;
849 triangle_badness*=priority;
850 /* mean triangle gradient */
851 triangle_gradient=ABS(Mz(tg->p0,M)-Mz(tg->p1,M))
852 + ABS(Mz(tg->p1,M)-Mz(tg->p2,M))
853 + ABS(Mz(tg->p2,M)-Mz(tg->p0,M));
854 /* multiply with priority value
855 * priority will change with the number of iterations */
856 /*triangle_gradient *=priority*0.3333; */
857
858
859 #ifdef DEBUG
860 printf("\n\nrefine_this: considering triangle %d points(%d,%d,%d) priority "MDF"\n",tg->n,tg->p0,tg->p1,tg->p2,priority);
861 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);
862 printf("refine_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges);
863 #endif
864 how_to_divide=0;
865 if ( fixed_points == 3 ) {
866 how_to_divide=3;
867 } else if ( override
868 || do_we_refine_this_triangle(area(tg),triangle_badness,triangle_gradient) ) {
869 /* find neighbour */
870 if ( (len2 >= len1) && (len2>=len0) ) {
871 tq.n=tg->t2;
872 neigh=RBT_find(&tq,&rt);
873
874 } else if ((len0 >= len1) && (len0>=len2) ) {
875 tq.n=tg->t0;
876 neigh=RBT_find(&tq,&rt);
877
878 } else {
879 tq.n=tg->t1;
880 neigh=RBT_find(&tq,&rt);
881 }
882 if ( (!neigh) || (neigh->status!=LIVE ) ) {
883 /* no neighbour exists, divide into 4 or 2 */
884 how_to_divide=4;
885 } else if ( neigh &&can_we_split_to_4(tg, neigh) ) {
886 /* if the largest edge is a boundary edge
887 * we know it cannot be flipped by delaunay refinement.
888 * so we split to 4
889 */
890 how_to_divide=4;
891 } else {
892 how_to_divide=3;
893 }
894 /* if the neighbours largest edge is also our largest edge,
895 * we can split to 4 */
896 }
897
898 if ( how_to_divide==3 ) {
899 #ifdef DEBUG
900 printf("subdivide_this: fixed_points=%d boundaries=%d\n",fixed_points,boundary_edges);
901 #endif
902 /* split to 3 by inserting point at centroid */
903 /* p0 /\ p0 /|\
904 * / \ / | \
905 * / \ to /1 | \ t1
906 * / \ t2 / / \3 \
907 * / \ / / pn\ \
908 * / \ / / \ \
909 * / p1 \ p2 // 2 \\
910 * -------------- ---------------
911 * p1 t0 p2
912 */
913 p.x=(Mx(tg->p0,M)+Mx(tg->p1,M)+Mx(tg->p2,M))*ONE_OVER_THREE ;
914 p.y=(My(tg->p0,M)+My(tg->p1,M)+My(tg->p2,M))*ONE_OVER_THREE;
915 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
916 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
917 exit(1);
918 }
919 p.z[0]=(Mz(tg->p0,M)+Mz(tg->p1,M)+Mz(tg->p2,M))*ONE_OVER_THREE;
920 p.val=VAR;
921
922 point_number=BIT_insert(&M,p);
923 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
924 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
925 exit(1);
926 }
927 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
928 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
929 exit(1);
930 }
931 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
932 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
933 exit(1);
934 }
935 tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
936 tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
937 tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
938 tg1->boundary=tg2->boundary=tg3->boundary=tg->boundary;
939 /* arrange neighbours */
940 tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
941 tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
942 tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
943 /* new triangles are alive */
944 tg1->status=tg2->status=tg3->status=LIVE;
945 /* old triangle is dead */
946 tg->status=DEAD;
947 #ifdef DEBUG
948 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);
949 #endif
950
951 /* update links from neighbours too */
952 tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
953 #ifdef DEBUG
954 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);
955 #endif
956 update_neighbours_around(tempt,tg->n,tg3->n);
957 tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
958 #ifdef DEBUG
959 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);
960 #endif
961 update_neighbours_around(tempt,tg->n,tg1->n);
962 tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
963 #ifdef DEBUG
964 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);
965 #endif
966 update_neighbours_around(tempt,tg->n,tg2->n);
967
968
969 parent_dag_node=tg->dag_p;
970 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
971 tg1->dag_p=current_dag_node;
972 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
973 tg2->dag_p=current_dag_node;
974 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
975 tg3->dag_p=current_dag_node;
976
977 RBT_insert((void*)tg1, &rt);
978 RBT_insert((void*)tg2, &rt);
979 RBT_insert((void*)tg3, &rt);
980 ntriangles+=3;
981
982 /* since all edges are on boundary, we cannot refine them */
983 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
984 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
985 delaunay_refine_constrained(tg3->n,point_number,&rt,&dt);
986
987 return(1);
988
989 } else if ( how_to_divide==4 ) {
990 #ifdef DEBUG
991 printf("subdivide_this: fixed_points=%d boundaries=%d badness="MDF" are="MDF"\n",fixed_points,boundary_edges,triangle_badness,area(tg));
992 #endif
993
994 /* split the largest edge by drawing a perpendicular */
995 /* general form */
996 /* p0 p0
997 * /\ /\
998 * / \ t3 / |\
999 * t0 / \ / | \
1000 * / t \ / 1| 2\
1001 * p1 /____p___\ p3 p1/____|___\ p3
1002 * \ / \ | 4 /
1003 * \ nei / \ 3| /
1004 * \ / \ | /
1005 * t1 \ / t2 \ |/
1006 * \/ \/
1007 * p2 p2
1008 */
1009 /* p1-p3 is the largest edge of triangle t */
1010 /* if p1-p3 is a boundary, neighbour is hidden and ignored */
1011 if ( (len2 >= len1) && (len2>=len0) ) {
1012 /* edge p0-p1 */
1013 p0=tg->p2;
1014 p1=tg->p0;
1015 p3=tg->p1;
1016 tq.n=tg->t2;
1017 neigh=RBT_find(&tq,&rt);
1018
1019 tq.n=tg->t1;
1020 t0=RBT_find(&tq,&rt);
1021 tq.n=tg->t0;
1022 t3=RBT_find(&tq,&rt);
1023 } else if ((len0 >= len1) && (len0>=len2) ) {
1024 /* edge p1-p2 */
1025 p0=tg->p0;
1026 p1=tg->p1;
1027 p3=tg->p2;
1028 tq.n=tg->t0;
1029 neigh=RBT_find(&tq,&rt);
1030
1031 tq.n=tg->t2;
1032 t0=RBT_find(&tq,&rt);
1033 tq.n=tg->t1;
1034 t3=RBT_find(&tq,&rt);
1035 } else {
1036 /* edge p2-p0 */
1037 p0=tg->p1;
1038 p1=tg->p2;
1039 p3=tg->p0;
1040 tq.n=tg->t1;
1041 neigh=RBT_find(&tq,&rt);
1042
1043 tq.n=tg->t0;
1044 t0=RBT_find(&tq,&rt);
1045 tq.n=tg->t2;
1046 t3=RBT_find(&tq,&rt);
1047 }
1048 /* determine foot of normal to longest edge */
1049 /* note the denominator will be zero only if the triangle is
1050 * a right angled one. but then the longest edge will be the diagonal
1051 * hence in practice it will never be zero unless points coincide */
1052 /* len0=-((Mx(p1)-Mx(p3))*(Mx(p3)-Mx(p0))+(My(p1)-My(p3))*(My(p3)-My(p0)))
1053 /((Mx(p1)-Mx(p0))*(Mx(p1)-Mx(p3))+(My(p1)-My(p0))*(My(p1)-My(p3)));
1054 p.x=(len0*Mx(p1)+Mx(p3))/(1+len0);
1055 p.y=(len0*My(p1)+My(p3))/(1+len0); */
1056 p.x=(Mx(p1,M)+Mx(p3,M))*0.5;
1057 p.y=(My(p1,M)+My(p3,M))*0.5;
1058 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
1059 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1060 exit(1);
1061 }
1062 /* p.z[0]=(len0*Mz(p1)+Mz(p3))/(1+len0); */
1063 p.z[0]=(Mz(p1,M)+Mz(p3,M))*0.5;
1064 point_number=is_on_an_edge(p1,p3);
1065 /* point_number= boundary_edge+1, if (p1,p3) lies on an edge */
1066 p.val=(point_number && edge_array[point_number-1].type==DR?FX:VAR);
1067 #ifdef DEBUG
1068 if (p.val==FX) {
1069 printf("subdivide_this: point fixed\n");
1070 } else {
1071 printf("subdivide_this: point variable\n");
1072 }
1073 #endif
1074 point_number=BIT_insert(&M,p);
1075 /* check for neighbour */
1076 if (neigh && neigh->status==LIVE) {
1077 if (neigh->t0==tg->n) {
1078 p2=neigh->p0;
1079 tq.n=neigh->t1;
1080 t1=RBT_find(&tq,&rt);
1081 tq.n=neigh->t2;
1082 t2=RBT_find(&tq,&rt);
1083 } else if (neigh->t1 == tg->n) {
1084 p2=neigh->p1;
1085 tq.n=neigh->t2;
1086 t1=RBT_find(&tq,&rt);
1087 tq.n=neigh->t0;
1088 t2=RBT_find(&tq,&rt);
1089 } else {
1090 p2=neigh->p2;
1091 tq.n=neigh->t0;
1092 t1=RBT_find(&tq,&rt);
1093 tq.n=neigh->t1;
1094 t2=RBT_find(&tq,&rt);
1095 }
1096 #ifdef DEBUG
1097 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);
1098 printf("sundivide_this: generalized to p0=%d,p1=%d,p2=%d,p3=%d\n",p0,p1,p2,p3);
1099 #endif
1100 }
1101 /*split to 4 or 2 */
1102 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
1103 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1104 exit(1);
1105 }
1106 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
1107 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1108 exit(1);
1109 }
1110 if ( neigh && neigh->status==LIVE) {
1111 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
1112 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1113 exit(1);
1114 }
1115 if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
1116 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1117 exit(1);
1118 }
1119 }
1120
1121 /* update triangle data */
1122 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
1123 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
1124 tg1->boundary=tg2->boundary=tg->boundary;
1125 /* alive, dead triangles */
1126 tg1->status=tg2->status=LIVE;
1127 tg->status=DEAD;
1128 if (neigh && neigh->status==LIVE) {
1129 tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
1130 tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;
1131 tg3->status=tg4->status=LIVE;
1132 tg3->boundary=tg4->boundary=neigh->boundary;
1133 /* at the last moment this */
1134 neigh->status=DEAD;
1135 tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
1136 tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
1137 tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
1138 tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
1139
1140 update_neighbours_around(t0,tg->n,tg1->n);
1141 update_neighbours_around(t3,tg->n,tg2->n);
1142 update_neighbours_around(t1,neigh->n,tg3->n);
1143 update_neighbours_around(t2,neigh->n,tg4->n);
1144 /* insert triangles into RBT */
1145 RBT_insert((void*)tg1,&rt);
1146 RBT_insert((void*)tg2,&rt);
1147 RBT_insert((void*)tg3,&rt);
1148 RBT_insert((void*)tg4,&rt);
1149 ntriangles+=4;
1150 /* update DAG */
1151 parent_dag_node=tg->dag_p;
1152 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1);
1153 tg1->dag_p=current_dag_node;
1154 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2);
1155 tg2->dag_p=current_dag_node;
1156 parent_dag_node=neigh->dag_p;
1157 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg3);
1158 tg3->dag_p=current_dag_node;
1159 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg4);
1160 tg4->dag_p=current_dag_node;
1161 #ifdef DEBUG
1162 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);
1163 #endif
1164 /* refine the triangulation */
1165 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
1166 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
1167 delaunay_refine_constrained(tg3->n,point_number,&rt,&dt);
1168 delaunay_refine_constrained(tg4->n,point_number,&rt,&dt);
1169
1170
1171 } else {
1172 tg1->t0=(neigh?neigh->n:-1);tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
1173 tg2->t0=(neigh?neigh->n:-1);tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
1174 update_neighbours_around(t0,tg->n,tg1->n);
1175 update_neighbours_around(t3,tg->n,tg2->n);
1176 /* insert triangles into RBT */
1177 RBT_insert((void*)tg1,&rt);
1178 RBT_insert((void*)tg2,&rt);
1179 ntriangles+=2;
1180 /* update DAG */
1181 parent_dag_node=tg->dag_p;
1182 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg1);
1183 tg1->dag_p=current_dag_node;
1184 current_dag_node=DAG_insert(&dt,parent_dag_node,(void*)tg2);
1185 tg2->dag_p=current_dag_node;
1186
1187 /* refine the triangulation */
1188 delaunay_refine_constrained(tg1->n,point_number,&rt,&dt);
1189 delaunay_refine_constrained(tg2->n,point_number,&rt,&dt);
1190
1191 }
1192 return(1);
1193
1194 }
1195
1196 return(0);
1197 }
1198
1199 /* similar to subdivide_all_triangles
1200 * but we take into account the gradient
1201 * at each point when dividing a triangle */
1202 int
refine_mesh_with_iterations(MY_DOUBLE priority)1203 refine_mesh_with_iterations(MY_DOUBLE priority)
1204 {
1205 triangle *rec;
1206 int flag=0;
1207
1208 /* gradient parameters */
1209 MY_DOUBLE grad_mean,max_grad,triangle_gradient;
1210 MY_INT count;
1211
1212 /* first do a traversal to see the max-min gradients are */
1213 grad_mean=max_grad=0;
1214 count=0;
1215 DAG_reverse_traverse_list_reset(&dt);
1216 rec=DAG_reverse_traverse_prune_list(&dt);
1217 while(rec){
1218 if (rec &&(rec->status==LIVE) ) {
1219 triangle_gradient=ABS(Mz(rec->p0,M)-Mz(rec->p1,M))
1220 + ABS(Mz(rec->p1,M)-Mz(rec->p2,M))
1221 + ABS(Mz(rec->p2,M)-Mz(rec->p0,M));
1222 if ( max_grad < triangle_gradient )
1223 max_grad=triangle_gradient;
1224 grad_mean=((MY_DOUBLE)count*grad_mean+triangle_gradient)/((MY_DOUBLE)count+1.0);
1225 count++;
1226 #ifdef DEBUG
1227 printf("refine_mesh_iter: %d gradient"MDF"\n",rec->n,triangle_gradient);
1228 #endif
1229 }
1230 rec=DAG_reverse_traverse_prune_list(&dt);
1231 }
1232
1233 /* refine all triangles with gradient >1/2(mean+max) */
1234 g_gradient_limit=0.5*(max_grad+grad_mean);
1235 #ifdef DEBUG
1236 printf("refine_mesh_iter: gradient max,mean "MDF"and "MDF"\n",max_grad,grad_mean);
1237 printf("refine_mesh_iter: set gradient limit"MDF"\n",g_gradient_limit);
1238 #endif
1239
1240 /* now do the proper refinement */
1241 DAG_reverse_traverse_list_reset(&dt);
1242 rec=DAG_reverse_traverse_prune_list(&dt);
1243 while(rec){
1244 if (rec &&(rec->status==LIVE)
1245 && refine_this_triangle(rec,0,priority)) {
1246 /* if a triangle was divided, flag=1 */
1247 #ifdef DEBUG
1248 printf("refining %d\n",rec->n);
1249 #endif
1250 flag |=1;
1251 }
1252 rec=DAG_reverse_traverse_prune_list(&dt);
1253 }
1254 return(flag);
1255 }
1256