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: pdnmesh.c,v 1.57 2005/04/24 05:32:24 sarod Exp $
17 */
18
19 #ifdef HAVE_CONFIG_H
20 #include <config.h>
21 #endif
22
23 #include <stdio.h>
24 #include <stdlib.h>
25 #include <string.h>
26 #include <ctype.h>
27 #include <math.h>
28 #ifndef WIN32
29 #include <unistd.h>
30 #include <gtk/gtk.h>
31 #include <gdk/gdkkeysyms.h>
32 #include <gtk/gtkgl.h>
33 #endif/*WIN32 */
34 #ifdef WIN32
35 #include <windows.h>
36 #endif/*!WIN32*/
37 #include <GL/gl.h>
38 #include <GL/glu.h>
39
40 #include "types.h"
41
42 int ntriangles,mouse_responce_flag,windowid;
43 DAG_tree dt;
44 RBT_tree rt;
45 double zoom_x1,zoom_x2,zoom_y1,zoom_y2;
46 GLdouble current_mat[16]; /* current zoom transform */
47 /* for getopt() */
48 extern char *optarg;
49 extern int optind, opterr, optopt;
50
51
52 int degree_of_freedom=1; /* degree of freedom for a point */
53 int requested_degree_of_freedom=0; /* user request, for eigenvalue problems */
54 int max_iteration_limit; /* limit for iterative solution */
55 int text_mode_only; /* flag to detect text mode operation */
56
57 /* boundary polygons */
58 extern boundary bound;
59
60 /* input cord filename */
61 char * input_cord_filename;
62
63 /* vector cross product */
64 #define CROSS_PRODUCT(x1,y1,x2,y2) \
65 ((x1)*(y2)-(y1)*(x2))
66 /* vector dot product */
67 #define DOT_PRODUCT(x1,y1,x2,y2) \
68 ((x1)*(x2)+(y1)*(y2))
69 /* function to test point p is on a line segment (p1-p2) */
70 static int
is_point_online_p(point p,point p1,point p2)71 is_point_online_p(point p, point p1,point p2)
72 {
73 /* p = p1 + lambda p2/(1+lambda) */
74 double lambda;
75
76 /*if (((p.x == p1.x)&&(p.y==p1.y))||((p.x==p2.x)&&(p.y==p2.y))) { */
77 if ((IS_ZERO(p.x - p1.x)&&IS_ZERO(p.y-p1.y))||(IS_ZERO(p.x-p2.x)&&IS_ZERO(p.y-p2.y))) {
78 /* points coincide==>on line */
79 return(1);
80 }
81 /* horizontal line */
82 if (IS_ZERO(p.x-p1.x)&&IS_ZERO(p.x-p2.x)) {
83 return((p.y-p1.y)/(p2.y-p.y) >0.0f);
84 }
85 /* vertical line */
86 if ( IS_ZERO(p.y-p1.y)&&IS_ZERO(p.y-p2.y) ) {
87 return((p.x-p1.x)/(p2.x-p.x) >0.0f);
88 }
89 lambda=(p.x-p1.x)/(p2.x-p.x);
90 if ( IS_ZERO(lambda - (p.y-p1.y)/(p2.y-p.y)) ) {
91 return(lambda>0.0f);
92 }
93 /* else */
94 return(0);
95 }
96
97 /* function to check point p is the same side as point c of line a-b */
98 /* return 1 if true, 0 if false, -1 if floating point over/under flow */
99 static int
same_side(point p,point c,point a,point b)100 same_side(point p, point c, point a, point b) {
101 double cp1,cp2;
102 cp1=CROSS_PRODUCT(b.x-a.x,b.y-a.y,p.x-a.x,p.y-a.y);
103 cp2=CROSS_PRODUCT(b.x-a.x,b.y-a.y,c.x-a.x,c.y-a.y);
104 #ifdef DEBUG
105 printf("[cross %lf,%lf",cp1,cp2);
106 #endif
107 /* else */
108 /*(cp1*cp2 >= 0) */
109 if ((cp1==0.0f)||(cp2==0.0f)
110 ||((cp1>0.0f)&&(cp2>0.0f))
111 ||((cp1<0.0f)&&(cp2<0.0f))) {
112 #ifdef DEBUG
113 printf("same side]");
114 #endif
115 return(1);
116 } else {
117 #ifdef DEBUG
118 printf("NOT same side]");
119 #endif
120 return(0);
121 }
122 }
123
124 static int
inclusion(const void * trig,const void * pt,const void * MM)125 inclusion(const void *trig,const void *pt, const void *MM/* Mesh*/){
126 triangle *t;
127 point *p,a,b,c;
128 Mesh *G;
129 t=(triangle *)trig;
130 p=(point *)pt;
131 G=(Mesh*)MM;
132 a.x=Mx(t->p0,*G);
133 a.y=My(t->p0,*G);
134 b.x=Mx(t->p1,*G);
135 b.y=My(t->p1,*G);
136 c.x=Mx(t->p2,*G);
137 c.y=My(t->p2,*G);
138 if (is_point_online_p(*p,a,b)
139 ||is_point_online_p(*p,b,c)
140 ||is_point_online_p(*p,c,a)) {
141 #ifdef DEBUG
142 printf("inclusion: point on boundary\n");
143 #endif
144 return(2);
145 }
146
147 /* if point is on line any triangle edge,
148 we consider that to be included */
149 #ifdef DEBUG
150 printf("inclusion (%lf,%lf)-(%lf,%lf)-(%lf,%lf) in (%lf,%lf)\n",
151 a.x,a.y,b.x,b.y,c.x,c.y,p->x,p->y);
152 #endif
153 /* check for errors */
154 if(same_side(*p,c,a,b)
155 &&same_side(*p,a,b,c)
156 &&same_side(*p,b,c,a)) {
157 /* we have floating point over/under flow
158 ans switch to alternative method */
159 #ifdef DEBUG
160 printf("inclusion: triangle %d includes the point\n",t->n);
161 #endif
162 return(1);
163 }
164 /* else */
165 #ifdef DEBUG
166 printf("inclusion: triangle %d NOT include the point\n",t->n);
167 #endif
168 return(0);
169 }
170
171
172 #ifdef EXTRA_SOURCE
173 static double
determ(double xx1,double yy1,double xx2,double yy2,double xx3,double yy3)174 determ( double xx1, double yy1, double xx2, double yy2, double xx3, double yy3 ) {
175 return ( xx1*(yy2-yy3)+xx2*(yy3-yy1)+xx3*(yy1-yy2) );
176 }
177
178 static int
ccw(point p0,point p1,point p2)179 ccw(point p0,point p1,point p2)
180 {
181 return(determ(p0.x,p0.y,p1.x,p1.y,p2.x,p2.y)>=0.0f);
182 }
183 /* is d inside triangle (a,b,c) ? */
184 #define INSIDE(a,b,c,d) \
185 (ccw(a,b,d)&&ccw(b,c,d)&&ccw(c,a,d))
186
187 /* checking inclusion of point in triangle */
188 static int
inclusion1(const void * trig,const void * pt,const void * MM)189 inclusion1(const void *trig,const void *pt, const void *MM)
190 {
191 triangle *t;
192 point *p,a,b,c;
193 Mesh *G;
194 t=(triangle *)trig;
195 p=(point *)pt;
196 G=(Mesh*)MM;
197 a.x=Mx(t->p0,*G);
198 a.y=My(t->p0,*G);
199 b.x=Mx(t->p1,*G);
200 b.y=My(t->p1,*G);
201 c.x=Mx(t->p2,*G);
202 c.y=My(t->p2,*G);
203 if (is_point_online_p(*p,a,b)
204 ||is_point_online_p(*p,b,c)
205 ||is_point_online_p(*p,c,a)) {
206 #ifdef DEBUG
207 printf("inclusion1: point on boundary\n");
208 #endif
209 return(2);
210 }
211 /* else */
212 return(INSIDE(a,b,c,*p));
213 }
214 #endif/* !EXTRA_SOURCE */
215
216 /**************************/
217 /* for triangles */
218 static int
equ(const void * a,const void * b)219 equ(const void*a,const void *b)
220 {
221 triangle *p,*q;
222 p=(triangle *)a;
223 q=(triangle *)b;
224
225 return(p->n==q->n);
226 }
227 static int
lt(const void * a,const void * b)228 lt(const void*a,const void *b)
229 {
230 triangle *p,*q;
231 p=(triangle *)a;
232 q=(triangle *)b;
233
234 return(p->n<q->n);
235
236 }
237
238 static void
print_detail(const void * a)239 print_detail(const void *a)
240 {
241 triangle *p;
242 p=(triangle *)a;
243 printf("%d :",p->n);
244 printf("%d:(%lf,%lf), %d:(%lf,%lf), %d:(%lf,%lf) ",p->p0,Mx(p->p0,M),My(p->p0,M),
245 p->p1,Mx(p->p1,M),My(p->p1,M),p->p2,Mx(p->p2,M),My(p->p2,M));
246 printf("neighbours %d,%d,%d ",p->t0,p->t1,p->t2);
247
248 p=(triangle *)p->dag_p->rec;
249 printf("parent in DAG points to %d\n",p->n);
250 }
251
252 static void
delete_triangle(void * a)253 delete_triangle(void *a)
254 {
255 triangle *p;
256 p=(triangle *)a;
257 if (!p) free(p);
258 }
259
260 /* a function to determine whether point p4 is inside */
261 /* the circumcircle of points p1,p2,p3 */
262 static int
inside_circle(int p1,int p2,int p3,int p4)263 inside_circle( int p1, int p2, int p3, int p4)
264 {
265 double temp1,temp2,temp3,temp4,a2,b2;
266 double A[3][3];
267 #ifdef DEBUG
268 printf("checking (%d,%d,%d) includes %d\n",p1,p2,p3,p4);
269 #endif
270 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));
271 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));
272 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));
273 /* we know this is not zero */
274 temp3 = 1/temp3;
275 a2 = (temp1*(My(p1,M)-My(p3,M))-temp2*(My(p1,M)-My(p2,M)))*temp3;
276 b2 = (temp2*(Mx(p1,M)-Mx(p2,M))-temp1*(Mx(p1,M)-Mx(p3,M)))*temp3;
277 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);
278 temp3=a2*(Mx(p1,M)-Mx(p4,M))+b2*(My(p1,M)-My(p4,M));
279 if ( ABS(temp3 - temp4) > TOL ) {
280 /* can make a simple decision */
281 if(temp3 < temp4 ) {
282 #ifdef DEBUG
283 printf("Simple case :points are inside circle\n");
284 #endif
285 return(1);
286 } else {
287 #ifdef DEBUG
288 printf("Simple case :points are NOT inside circle\n");
289 #endif
290 return(0);
291 }
292 } else {
293 /* this is a close call. evaluate determinant */
294 A[0][0]=Mx(p1,M)-Mx(p4,M);A[0][1]=My(p1,M)-My(p4,M);A[0][2]=A[0][0]*A[0][0]+A[0][1]*A[0][1];
295 A[1][0]=Mx(p2,M)-Mx(p4,M);A[1][1]=My(p2,M)-My(p4,M);A[1][2]=A[1][0]*A[1][0]+A[1][1]*A[1][1];
296 A[2][0]=Mx(p3,M)-Mx(p4,M);A[2][1]=My(p3,M)-My(p4,M);A[2][2]=A[2][0]*A[2][0]+A[2][1]*A[2][1];
297 /* evaluate determinant */
298 temp3=A[0][0]*(A[1][2]*A[2][2]-A[2][1]*A[1][2])
299 -A[0][1]*(A[1][0]*A[2][2]-A[2][0]*A[1][2])
300 +A[0][2]*(A[1][0]*A[2][1]-A[2][0]*A[1][1]);
301 if(temp3 <=0) {
302 #ifdef DEBUG
303 printf("Hard case :points are inside circle\n");
304 #endif
305 return(1);
306 } else {
307 #ifdef DEBUG
308 printf("Hard case :points are NOT inside circle\n");
309 #endif
310 return(0);
311 }
312 }
313 return(0);
314 }
315
316 /* update neighbour link in a triangle */
317 /* neighbour old will be replaced by neighbour new */
318 static void
update_neighbours_around(triangle * neighbour,int old,int new)319 update_neighbours_around(triangle *neighbour,int old,int new)
320 {
321 if ( neighbour ) {
322 if ( neighbour->t0 == old ) {
323 neighbour->t0 = new;
324 } else if ( neighbour->t1 == old ) {
325 neighbour->t1 = new;
326 } else if ( neighbour->t2 == old ) {
327 neighbour->t2 = new;
328 }
329 }
330 }
331
332 /* tgn - triangle number to refine */
333 /* ptn - point number of vertex of tgn facing the edge to be refined */
334 /* assertion - at each step, the remainder of triangles form a valid
335 * delaunay triangulation */
336 int
delaunay_refine(int tgn,int ptn,RBT_tree * rbt,DAG_tree * dag)337 delaunay_refine(int tgn, int ptn, RBT_tree * rbt,DAG_tree *dag)
338 {
339 triangle tq,*t,*t_n;
340 int neighbour,p0,p1,p2,p4;
341 int t1,t2,t3,t4; /* neighbours */
342 triangle *tp1,*tp2,*tp3,*tp4,*tnew1,*tnew2;
343
344 DAG_node *current_dag_node,*parent_dag_node;
345 /* fix everything into this generalized form
346 * p2/\
347 * /| \
348 * / | \
349 * t1/ | \ t4
350 * / | \
351 * p0 / | \ p4
352 * \ t | tn /
353 * \ | /
354 * t2 \ | / t3
355 * \ | /
356 * p1
357 */
358
359 /* get the triangle */
360 tq.n=tgn;
361 t=RBT_find(&tq, rbt);
362 if ( t && ( t->n == tgn )) { /* triangle exists */
363 /* first find neighbour facing the edge to be refined */
364 if ( ptn == t->p0 ) {
365 neighbour=t->t0;
366 p0=t->p0;p1=t->p1;p2=t->p2;
367 t1=t->t1; t2=t->t2;
368 } else if ( ptn == t->p1 ) {
369 neighbour=t->t1;
370 p0=t->p1;p1=t->p2;p2=t->p0;
371 t1=t->t2; t2=t->t0;
372 } else {
373 neighbour=t->t2;
374 p0=t->p2;p1=t->p0;p2=t->p1;
375 t1=t->t0; t2=t->t1;
376 }
377 /* if the edge to be flipped is a fixed one,
378 * we cannot do anything, we quit */
379 /*if ((Mval(p1)==FX) &&(Mval(p2)==FX)){
380 printf("cannot refine\n");
381 return(1);
382 }*/
383 /* now find the neighbouring triangle */
384 /* it may not exist */
385 tq.n=neighbour;
386 t_n=RBT_find(&tq, rbt);
387 if ( !t_n ) {/* no neighbour */
388 return(0);
389 }
390 /* find the shared edge with neigbour, or rather, the opposing vertex */
391 if ( tgn == t_n->t0 ) {
392 p4=t_n->p0;
393 t3=t_n->t1; t4=t_n->t2;
394 } else if ( tgn == t_n->t1 ) {
395 p4=t_n->p1;
396 t3=t_n->t2; t4=t_n->t0;
397 } else {
398 p4=t_n->p2;
399 t3=t_n->t0; t4=t_n->t1;
400 }
401 /* now we have fitted everything to the general framework */
402 if ( inside_circle(p0,p2,p1,p4) ) {
403 /* we flip the edge (p1,p2) */
404 /* this could only be done if p4 can be seen
405 * by p0, i.e. (p1,p2) is the largest edge of neighbour
406 * but since we assert a valid initial triangulation,
407 * this will hold.
408 */
409
410 /* get neighbouring triangles first */
411 tq.n=t1;tp1=RBT_find(&tq, rbt);
412 tq.n=t2;tp2=RBT_find(&tq, rbt);
413 tq.n=t3;tp3=RBT_find(&tq, rbt);
414 tq.n=t4;tp4=RBT_find(&tq, rbt);
415 /* create 2 new triangles */
416 if ((tnew1=(triangle *)malloc(sizeof(triangle)))==0) {
417 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
418 exit(1);
419 }
420 if ((tnew2=(triangle *)malloc(sizeof(triangle)))==0) {
421 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
422 exit(1);
423 }
424 /* we transform the setup into
425 * p2/\
426 * / \
427 * / \
428 * t1/ \ t4
429 * / new1 \
430 * p0 / \ p4
431 * \- - ------ /
432 * \ new2 /
433 * t2 \ / t3
434 * \ /
435 * p1
436 */
437 tnew1->p0=p0;tnew1->p1=p4;tnew1->p2=p2;
438 tnew1->n=ntriangles;
439 tnew2->p0=p0;tnew2->p1=p1;tnew2->p2=p4;
440 tnew2->n=ntriangles+1;
441 tnew1->status=tnew2->status=LIVE;
442 t->status=t_n->status=DEAD;
443 /* neighbours */
444 tnew1->t0=t4;tnew1->t1=t1;tnew1->t2=tnew2->n;
445 tnew2->t0=t3;tnew2->t1=tnew1->n;tnew2->t2=t2;
446 /* undate the neighbours of the neighbours too */
447 update_neighbours_around(tp1,tgn,tnew1->n);
448 update_neighbours_around(tp2,tgn,tnew2->n);
449 update_neighbours_around(tp3,neighbour,tnew2->n);
450 update_neighbours_around(tp4,neighbour,tnew1->n);
451 /* insert the 2 new triangles into RB tree */
452 RBT_insert((void*)tnew1, rbt);
453
454 parent_dag_node=t->dag_p;
455 current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew1);
456 tnew1->dag_p=current_dag_node;
457 /* link it to the other parent in dag */
458 DAG_link(dag,t_n->dag_p,current_dag_node);
459
460 /* second triangle */
461 RBT_insert((void*)tnew2, rbt);
462
463 parent_dag_node=t->dag_p;
464 current_dag_node=DAG_insert(dag, parent_dag_node, (void *)tnew2);
465 tnew2->dag_p=current_dag_node;
466 /* link it to the other parent in dag */
467 DAG_link(dag,t_n->dag_p,current_dag_node);
468
469 /* increment triangle count */
470 ntriangles +=2;
471 /* printf("flipped %d neighbour %d to %d and %d\n",tgn,neighbour,tnew1->n,tnew2->n); */
472
473 /* call this function recursively */
474 delaunay_refine(tnew1->n,p0,rbt,dag);
475 delaunay_refine(tnew2->n,p0,rbt,dag);
476 }
477
478
479 } else {
480 /* error */
481 printf("cannot find triangle\n");
482 return(-1);
483 }
484
485
486 return(0);
487 }
488
489
490 /* function to test point is on a line segment */
491 static double
is_point_online(int p,int p1,int p2)492 is_point_online(int p, int p1,int p2)
493 {
494 /* p = p1 + lambda p2/(1+lambda) */
495 double x,x1,x2,y,y1,y2,lambda;
496
497 x=Mx(p,M);
498 x1=Mx(p1,M);
499 x2=Mx(p2,M);
500 y=My(p,M);
501 y1=My(p1,M);
502 y2=My(p2,M);
503
504 /*printf("(%lf,%lf)(%lf,%lf)(%lf,%lf)\n",x,y,x1,y1,x2,y2); */
505 if ((IS_ZERO(x - x1)&&IS_ZERO(y-y1))||(IS_ZERO(x-x2)&&IS_ZERO(y-y2))) {
506 /* points coincide */
507 return(-1);
508 }
509 if (IS_ZERO(x-x1)&&IS_ZERO(x-x2)) {
510 return((y-y1)/(y2-y));
511 }
512 if ( IS_ZERO(y-y1)&&IS_ZERO(y-y2) ) {
513 return((x-x1)/(x2-x));
514 }
515 lambda=(x-x1)/(x2-x);
516 if ( IS_ZERO(lambda - (y-y1)/(y2-y)) ) {
517 return(lambda);
518 }
519 /*printf("l1=%lf,l2=%lf\n",(x-x1)/(x2-x),(y-y1)/(y2-y)); */
520 /* else */
521 return(-1);
522 }
523
524 int
insert_point_to_mesh(point p)525 insert_point_to_mesh(point p)
526 {
527 DAG_node *current_dag_node,*parent_dag_node;
528 triangle *tg,*tg1,*tg2,*tg3,*tempt,tq;
529 triangle *tg4,*t0,*t1,*t2,*t3,*neigh;
530 int status,p0,p1,p2,p3,point_number,neighbour_number;
531 double lambda1,lambda2,lambda3;
532
533
534 point_number=BIT_insert(&M,p);
535 #ifdef DEBUG
536 printf("insert_point_to_mesh: new point %d\n",point_number);
537 printf("insert_point_to_mesh: Insert point count=%d, maximum=%d\n",M.count,M.maxpoints);
538 #endif
539 current_dag_node=DAG_find(&dt, (void *)&p,(void*)&M);
540 if ( current_dag_node ) {
541 /* point is included in this triangle */
542 tg=(triangle *)current_dag_node->rec;
543 #ifdef DEBUG
544 printf("insert_p_..:point inside %d->(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2);
545 #endif
546
547 /* insert point into bitree */
548 lambda1=is_point_online(point_number,tg->p0,tg->p1);
549 lambda2=is_point_online(point_number,tg->p1,tg->p2);
550 lambda3=is_point_online(point_number,tg->p2,tg->p0);
551 #ifdef DEBUG
552 printf("insert_p_..:lambdas %lf,%lf,%lf\n",lambda1,lambda2,lambda3);
553 #endif
554 if ( (lambda1<0)&&(lambda2<0)&&(lambda3<0)) {
555
556 /* split the current triangle into 3 */
557 /* p0 /\ p0 /|\
558 * / \ / | \
559 * / \ to /1 | \ t1
560 * / \ t2 / / \3 \
561 * / \ / / pn\ \
562 * / \ / / \ \
563 * / p1 \ p2 // 2 \\
564 * -------------- ---------------
565 * p1 t0 p2
566 */
567 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
568 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
569 exit(1);
570 }
571 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
572 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
573 exit(1);
574 }
575 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
576 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
577 exit(1);
578 }
579 tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
580 tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
581 tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
582 /* arrange neighbours */
583 tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
584 tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
585 tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
586 /* new triangles are alive */
587 tg1->status=tg2->status=tg3->status=LIVE;
588 /* old triangle is dead */
589 tg->status=DEAD;
590
591 /* update links from neighbours too */
592 tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
593 update_neighbours_around(tempt,tg->n,tg3->n);
594 tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
595 update_neighbours_around(tempt,tg->n,tg1->n);
596 tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
597 update_neighbours_around(tempt,tg->n,tg2->n);
598
599
600 parent_dag_node=current_dag_node;
601 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
602 tg1->dag_p=current_dag_node;
603 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
604 tg2->dag_p=current_dag_node;
605 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
606 tg3->dag_p=current_dag_node;
607
608 status = RBT_insert((void*)tg1, &rt);
609 status = RBT_insert((void*)tg2, &rt);
610 status = RBT_insert((void*)tg3, &rt);
611 ntriangles+=3;
612
613 /* now Delaunay refine the neighbouring triangles */
614 delaunay_refine(tg1->n,point_number,&rt,&dt);
615 delaunay_refine(tg2->n,point_number,&rt,&dt);
616 delaunay_refine(tg3->n,point_number,&rt,&dt);
617 } else {
618 /* one lambda is positive */
619 /* split to 4 */
620 /* general form */
621 /* p0 p0
622 * /\ /\
623 * / \ t3 / |\
624 * t0 / \ / | \
625 * / t \ / 1| 2\
626 * p1 /____p___\ p3 p1/____|___\ p3
627 * \ / \ | 4 /
628 * \ nei / \ 3| /
629 * \ / \ | /
630 * t1 \ / t2 \ |/
631 * \/ \/
632 * p2 p2
633 */
634 /*printf("split 4\n"); */
635 neigh=0;
636 if ( lambda1 > 0 ) {
637 p0=tg->p2;
638 p1=tg->p0;
639 p3=tg->p1;
640 tq.n=tg->t2;
641 neighbour_number=tq.n;
642 neigh=RBT_find(&tq, &rt);
643
644 tq.n=tg->t1;
645 t0=RBT_find(&tq, &rt);
646 tq.n=tg->t0;
647 t3=RBT_find(&tq, &rt);
648
649 } else if ( lambda2 > 0 ) {
650 p0=tg->p0;
651 p1=tg->p1;
652 p3=tg->p2;
653 tq.n=tg->t0;
654 neighbour_number=tq.n;
655 neigh=RBT_find(&tq, &rt);
656
657 tq.n=tg->t2;
658 t0=RBT_find(&tq, &rt);
659 tq.n=tg->t1;
660 t3=RBT_find(&tq, &rt);
661
662 } else { /* lambda3 > 0 */
663 p0=tg->p1;
664 p1=tg->p2;
665 p3=tg->p0;
666 tq.n=tg->t1;
667 neighbour_number=tq.n;
668 neigh=RBT_find(&tq, &rt);
669
670 tq.n=tg->t0;
671 t0=RBT_find(&tq, &rt);
672 tq.n=tg->t2;
673 t3=RBT_find(&tq, &rt);
674
675 }
676 /* if neighbour exists */
677 if ( neigh ) {
678 /*printf("neighbour exists %d\n",neigh->n); */
679 if ( neigh->t0 == tg->n ) {
680 p2=neigh->p0;
681 tq.n=neigh->t1;
682 t1=RBT_find(&tq,&rt);
683 tq.n=neigh->t2;
684 t2=RBT_find(&tq,&rt);
685 } else if (neigh->t1 == tg->n ) {
686 p2=neigh->p1;
687 tq.n=neigh->t2;
688 t1=RBT_find(&tq,&rt);
689 tq.n=neigh->t0;
690 t2=RBT_find(&tq,&rt);
691 } else {
692 p2=neigh->p2;
693 tq.n=neigh->t0;
694 t1=RBT_find(&tq,&rt);
695 tq.n=neigh->t1;
696 t2=RBT_find(&tq,&rt);
697 }
698
699 /* printf("generalized %d %d %d %d\n",p0,p1,p2,p3); */
700 /* split to 4 */
701 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
702 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
703 exit(1);
704 }
705 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
706 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
707 exit(1);
708 }
709 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
710 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
711 exit(1);
712 }
713 if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
714 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
715 exit(1);
716 }
717 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
718 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
719 tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
720 tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;
721
722 /* alive, dead triangles */
723 /* new triangles are alive */
724 tg1->status=tg2->status=tg3->status=tg4->status=LIVE;
725 /* old triangle is dead */
726 tg->status=neigh->status=DEAD;
727
728 /* arrange neighbours */
729 tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0 ? t0->n:-1);
730 tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
731 tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
732 tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
733 /* update old neighbours */
734 update_neighbours_around(t0,tg->n,tg1->n);
735 update_neighbours_around(t3,tg->n,tg2->n);
736 update_neighbours_around(t1,neigh->n,tg3->n);
737 update_neighbours_around(t2,neigh->n,tg4->n);
738
739 /* insert into RBT */
740 RBT_insert((void*)tg1, &rt);
741 RBT_insert((void*)tg2, &rt);
742 RBT_insert((void*)tg3, &rt);
743 RBT_insert((void*)tg4, &rt);
744 ntriangles+=4;
745
746 /* update DAG */
747 parent_dag_node=tg->dag_p;
748 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
749 tg1->dag_p=current_dag_node;
750 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
751 tg2->dag_p=current_dag_node;
752 parent_dag_node=neigh->dag_p;
753 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
754 tg3->dag_p=current_dag_node;
755 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg4);
756 tg4->dag_p=current_dag_node;
757
758 delaunay_refine(tg1->n,point_number,&rt,&dt);
759 delaunay_refine(tg2->n,point_number,&rt,&dt);
760 delaunay_refine(tg3->n,point_number,&rt,&dt);
761 delaunay_refine(tg4->n,point_number,&rt,&dt);
762
763 /* printf("split to %d %d %d %d\n",tg1->n,tg2->n,tg3->n,tg4->n);
764 DAG_traverse(&dt); */
765 } else { /* no neighbour */
766 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
767 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
768 exit(1);
769 }
770 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
771 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
772 exit(1);
773 }
774
775 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
776 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
777
778 /* alive, dead triangles */
779 /* new triangles are alive */
780 tg1->status=tg2->status=LIVE;
781 /* old triangle is dead */
782 tg->status=DEAD;
783
784 /* arrange neighbours */
785 tg1->t0=neighbour_number;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
786 tg2->t0=neighbour_number;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
787 /* update old neighbours */
788 update_neighbours_around(t0,tg->n,tg1->n);
789 update_neighbours_around(t3,tg->n,tg2->n);
790
791 /* insert into RBT */
792 RBT_insert((void*)tg1, &rt);
793 RBT_insert((void*)tg2, &rt);
794 ntriangles+=2;
795
796 /* update DAG */
797 parent_dag_node=tg->dag_p;
798 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
799 tg1->dag_p=current_dag_node;
800 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
801 tg2->dag_p=current_dag_node;
802 delaunay_refine(tg1->n,point_number,&rt,&dt);
803 delaunay_refine(tg2->n,point_number,&rt,&dt);
804
805 }
806 }
807
808 }
809 return(point_number);
810 }
811
812 /***************************/
813 /* remove (not literally) triangles outside a polygon */
814 static void
remove_exterior_and_update_triangles(boundary * b)815 remove_exterior_and_update_triangles(boundary *b)
816 {
817 triangle *rec;
818 point p;
819 int i;
820 #ifdef DEBUG
821 int j;
822 triangle brec;
823 #endif
824 /* remove outer boundary triangles */
825 DAG_traverse_list_reset(&dt);
826 rec=DAG_traverse_prune_list(&dt);
827 while (rec) {
828 if (rec &&rec->status==LIVE) {
829 p.x=(Mx(rec->p0,M)+Mx(rec->p1,M)+Mx(rec->p2,M))/3;
830 p.y=(My(rec->p0,M)+My(rec->p1,M)+My(rec->p2,M))/3;
831 if (!inside_poly_point(&(b->poly_array[0]),p)) {
832 rec->status=HIDDEN;
833 #ifdef DEBUG
834 printf("remove_ext: updated triangle %d(%d,%d,%d) status %d\n",rec->n,rec->p0,rec->p1,rec->p2,rec->status);
835 #endif
836
837
838 }
839 }
840
841 rec=DAG_traverse_prune_list(&dt);
842 }
843
844 /* update other boundary triangles */
845 for (i=0;i<b->npoly;i++) {
846 DAG_traverse_list_reset(&dt);
847 rec=DAG_traverse_prune_list(&dt);
848 while (rec) {
849 if (rec &&rec->status==LIVE) {
850 p.x=(Mx(rec->p0,M)+Mx(rec->p1,M)+Mx(rec->p2,M))/3;
851 p.y=(My(rec->p0,M)+My(rec->p1,M)+My(rec->p2,M))/3;
852 if (inside_poly_point(&(b->poly_array[i]),p)) {
853 rec->boundary=i;
854 /* if hollow boundary, make triangle hidden*/
855 if ( (b->poly_array[i]).hollow==0 ) {
856 rec->status=HIDDEN;
857 }
858 }
859 }
860 rec=DAG_traverse_prune_list(&dt);
861 }
862 }
863
864 #ifdef DEBUG
865 printf("===============\n");
866 for (j=0;j<ntriangles;j++) {
867 brec.n=j;
868 rec= RBT_find(&brec, &rt);
869
870 if (rec->status!=DEAD) {
871 printf("remove_ext: check triangle %d(%d,%d,%d) status %d\n",rec->n,rec->p0,rec->p1,rec->p2,rec->status);
872 }
873 }
874 #endif
875 }
876
877 /*************************/
878 /* generate mesh */
879 void
generate_mesh(const char * filename)880 generate_mesh(const char *filename) {
881
882 DAG_node *current_dag_node,*parent_dag_node;
883 triangle *tg,*tg1,*tg2,*tg3,*tempt,tq;
884 triangle *tg4,*t0,*t1,*t2,*t3,*neigh;
885 int status,p0,p1,p2,p3,total,i,point_number;
886 point p;
887 point *point_array;
888 double lambda1,lambda2,lambda3;
889 polygon temp_polygon;
890 poly_edge e;
891
892 /* read input */
893 total=read_input_file(&point_array,filename);
894
895 if ((tg=(triangle *)malloc(sizeof(triangle)))==0){
896 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
897 exit(1);
898 }
899 p.x=point_array[0].x;p.y=point_array[0].y;
900 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
901 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
902 exit(1);
903 }
904 p.z[0]=0;
905 p.val=FX; /* we dont care they are FX or not */
906 #ifdef DEBUG
907 printf("inserting point %d "MDF","MDF"\n",0,p.x,p.y);
908 #endif
909 tg->p0=BIT_insert(&M, p );
910 p.x=point_array[1].x;p.y=point_array[1].y;
911 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
912 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
913 exit(1);
914 }
915 p.z[0]=0;
916 p.val=FX;
917 #ifdef DEBUG
918 printf("inserting point %d "MDF","MDF"\n",1,p.x,p.y);
919 #endif
920 tg->p1=BIT_insert(&M, p );
921 p.x=point_array[2].x;p.y=point_array[2].y;
922 if ((p.z=(MY_DOUBLE*)malloc((size_t)(degree_of_freedom)*sizeof(MY_DOUBLE)))==0) {
923 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
924 exit(1);
925 }
926 p.z[0]=0;
927 p.val=FX;
928 #ifdef DEBUG
929 printf("inserting point %d "MDF","MDF"\n",2,p.x,p.y);
930 #endif
931 tg->p2=BIT_insert(&M, p );
932
933 tg->n=0;
934 /* triangle is alive */
935 tg->status=LIVE;
936 /* make neighbours -1 */
937 tg->t0=tg->t1=tg->t2=-1;
938 p0=tg->p0;p1=tg->p1;p2=tg->p2;
939 ntriangles=1;
940
941 status = RBT_insert((void*)tg, &rt);
942 parent_dag_node=dt.root;
943 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg);
944 /* set pointer to dag */
945 tg->dag_p=current_dag_node;
946 /* we have inserted the outer triangle */
947 /* now process the remaining points */
948 for (i=3;i<total;i++) {
949 p.x=point_array[i].x;p.y=point_array[i].y;p.val=point_array[i].val;
950 p.z=point_array[i].z;
951 point_number=BIT_insert(&M,p);
952 #ifdef DEBUG
953 printf("do_intial_triangulation: inserting point %d "MDF","MDF" z="MDF" as %d\n",i,p.x,p.y,*p.z,point_number);
954 #endif
955 if ( point_number != i) {
956 fprintf(stderr,"***: insertion of new point %d returned an old point %d.\n your input file is screwed up. you may have duplicate points. aborting..\n",i,point_number);
957 abort();
958 }
959 current_dag_node=DAG_find(&dt, (void *)&p,(void*)&M);
960 if ( current_dag_node ) {
961 /* point is included in this triangle */
962 tg=(triangle *)current_dag_node->rec;
963 #ifdef DEBUG
964 printf("do_intial_triangulation: point inside %d->(%d,%d,%d)\n",tg->n,tg->p0,tg->p1,tg->p2);
965 #endif
966 /* insert point into bitree */
967 lambda1=is_point_online(point_number,tg->p0,tg->p1);
968 lambda2=is_point_online(point_number,tg->p1,tg->p2);
969 lambda3=is_point_online(point_number,tg->p2,tg->p0);
970 #ifdef DEBUG
971 printf("lambdas %lf,%lf,%lf\n",lambda1,lambda2,lambda3);
972 #endif
973 if ( (lambda1<0)&&(lambda2<0)&&(lambda3<0)) {
974
975 /* split the current triangle into 3 */
976 /* p0 /\ p0 /|\
977 * / \ / | \
978 * / \ to /1 | \ t1
979 * / \ t2 / / \3 \
980 * / \ / / pn\ \
981 * / \ / / \ \
982 * / p1 \ p2 // 2 \\
983 * -------------- ---------------
984 * p1 t0 p2
985 */
986 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
987 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
988 exit(1);
989 }
990 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
991 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
992 exit(1);
993 }
994 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
995 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
996 exit(1);
997 }
998 tg1->p0=tg->p0;tg1->p1=tg->p1;tg1->p2=point_number;tg1->n=ntriangles;
999 tg2->p0=tg->p1;tg2->p1=tg->p2;tg2->p2=point_number;tg2->n=ntriangles+1;
1000 tg3->p0=tg->p2;tg3->p1=tg->p0;tg3->p2=point_number;tg3->n=ntriangles+2;
1001 /* arrange neighbours */
1002 tg1->t0=tg2->n;tg1->t1=tg3->n;tg1->t2=tg->t2;
1003 tg2->t0=tg3->n;tg2->t1=tg1->n;tg2->t2=tg->t0;
1004 tg3->t0=tg1->n;tg3->t1=tg2->n;tg3->t2=tg->t1;
1005 /* new triangles are alive */
1006 tg1->status=tg2->status=tg3->status=LIVE;
1007 /* old triangle is dead */
1008 tg->status=DEAD;
1009
1010 /* update links from neighbours too */
1011 tq.n=tg->t1;tempt=RBT_find(&tq, &rt);
1012 update_neighbours_around(tempt,tg->n,tg3->n);
1013 tq.n=tg->t2;tempt=RBT_find(&tq, &rt);
1014 update_neighbours_around(tempt,tg->n,tg1->n);
1015 tq.n=tg->t0;tempt=RBT_find(&tq, &rt);
1016 update_neighbours_around(tempt,tg->n,tg2->n);
1017
1018
1019 parent_dag_node=current_dag_node;
1020 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
1021 tg1->dag_p=current_dag_node;
1022 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
1023 tg2->dag_p=current_dag_node;
1024 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
1025 tg3->dag_p=current_dag_node;
1026
1027 status = RBT_insert((void*)tg1, &rt);
1028 status = RBT_insert((void*)tg2, &rt);
1029 status = RBT_insert((void*)tg3, &rt);
1030 ntriangles+=3;
1031
1032 /* now Delaunay refine the neighbouring triangles */
1033 delaunay_refine(tg1->n,point_number,&rt,&dt);
1034 delaunay_refine(tg2->n,point_number,&rt,&dt);
1035 delaunay_refine(tg3->n,point_number,&rt,&dt);
1036 } else {
1037 /* one lambda is positive */
1038 /* split to 4 */
1039 /* general form */
1040 /* p0 p0
1041 * /\ /\
1042 * / \ t3 / |\
1043 * t0 / \ / | \
1044 * / t \ / 1| 2\
1045 * p1 /____p___\ p3 p1/____|___\ p3
1046 * \ / \ | 4 /
1047 * \ nei / \ 3| /
1048 * \ / \ | /
1049 * t1 \ / t2 \ |/
1050 * \/ \/
1051 * p2 p2
1052 */
1053 /*printf("split 4\n"); */
1054 neigh=0;
1055 if ( lambda1 > 0 ) {
1056 p0=tg->p2;
1057 p1=tg->p0;
1058 p3=tg->p1;
1059 tq.n=tg->t2;
1060 neigh=RBT_find(&tq, &rt);
1061
1062 tq.n=tg->t1;
1063 t0=RBT_find(&tq, &rt);
1064 tq.n=tg->t0;
1065 t3=RBT_find(&tq, &rt);
1066
1067 } else if ( lambda2 > 0 ) {
1068 p0=tg->p0;
1069 p1=tg->p1;
1070 p3=tg->p2;
1071 tq.n=tg->t0;
1072 neigh=RBT_find(&tq, &rt);
1073
1074 tq.n=tg->t2;
1075 t0=RBT_find(&tq, &rt);
1076 tq.n=tg->t1;
1077 t3=RBT_find(&tq, &rt);
1078
1079 } else { /* lambda3 > 0 */
1080 p0=tg->p1;
1081 p1=tg->p2;
1082 p3=tg->p0;
1083 tq.n=tg->t1;
1084 neigh=RBT_find(&tq, &rt);
1085
1086 tq.n=tg->t0;
1087 t0=RBT_find(&tq, &rt);
1088 tq.n=tg->t2;
1089 t3=RBT_find(&tq, &rt);
1090
1091 }
1092 /* if neighbour exists */
1093 if ( neigh ) {
1094 /*printf("neighbour exists %d\n",neigh->n); */
1095 if ( neigh->t0 == tg->n ) {
1096 p2=neigh->p0;
1097 tq.n=neigh->t1;
1098 t1=RBT_find(&tq,&rt);
1099 tq.n=neigh->t2;
1100 t2=RBT_find(&tq,&rt);
1101 } else if (neigh->t1 == tg->n ) {
1102 p2=neigh->p1;
1103 tq.n=neigh->t2;
1104 t1=RBT_find(&tq,&rt);
1105 tq.n=neigh->t0;
1106 t2=RBT_find(&tq,&rt);
1107 } else {
1108 p2=neigh->p2;
1109 tq.n=neigh->t0;
1110 t1=RBT_find(&tq,&rt);
1111 tq.n=neigh->t1;
1112 t2=RBT_find(&tq,&rt);
1113 }
1114
1115 /* printf("generalized %d %d %d %d\n",p0,p1,p2,p3); */
1116 /* split to 4 */
1117 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
1118 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1119 exit(1);
1120 }
1121 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
1122 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1123 exit(1);
1124 }
1125 if ((tg3=(triangle *)malloc(sizeof(triangle)))==0) {
1126 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1127 exit(1);
1128 }
1129 if ((tg4=(triangle *)malloc(sizeof(triangle)))==0) {
1130 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1131 exit(1);
1132 }
1133 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
1134 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
1135 tg3->p0=point_number;tg3->p1=p1;tg3->p2=p2;tg3->n=ntriangles+2;
1136 tg4->p0=point_number;tg4->p1=p2;tg4->p2=p3;tg4->n=ntriangles+3;
1137
1138 /* alive, dead triangles */
1139 /* new triangles are alive */
1140 tg1->status=tg2->status=tg3->status=tg4->status=LIVE;
1141 /* old triangle is dead */
1142 tg->status=neigh->status=DEAD;
1143
1144 /* arrange neighbours */
1145 tg1->t0=tg3->n;tg1->t1=tg2->n;tg1->t2=(t0 ? t0->n:-1);
1146 tg2->t0=tg4->n;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
1147 tg3->t0=(t1?t1->n:-1);tg3->t1=tg4->n;tg3->t2=tg1->n;
1148 tg4->t0=(t2?t2->n:-1);tg4->t1=tg2->n;tg4->t2=tg3->n;
1149 /* update old neighbours */
1150 update_neighbours_around(t0,tg->n,tg1->n);
1151 update_neighbours_around(t3,tg->n,tg2->n);
1152 update_neighbours_around(t1,neigh->n,tg3->n);
1153 update_neighbours_around(t2,neigh->n,tg4->n);
1154
1155 /* insert into RBT */
1156 RBT_insert((void*)tg1, &rt);
1157 RBT_insert((void*)tg2, &rt);
1158 RBT_insert((void*)tg3, &rt);
1159 RBT_insert((void*)tg4, &rt);
1160 ntriangles+=4;
1161
1162 /* update DAG */
1163 parent_dag_node=tg->dag_p;
1164 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
1165 tg1->dag_p=current_dag_node;
1166 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
1167 tg2->dag_p=current_dag_node;
1168 parent_dag_node=neigh->dag_p;
1169 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg3);
1170 tg3->dag_p=current_dag_node;
1171 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg4);
1172 tg4->dag_p=current_dag_node;
1173
1174 delaunay_refine(tg1->n,point_number,&rt,&dt);
1175 delaunay_refine(tg2->n,point_number,&rt,&dt);
1176 delaunay_refine(tg3->n,point_number,&rt,&dt);
1177 delaunay_refine(tg4->n,point_number,&rt,&dt);
1178
1179 /*printf("split to %d %d %d %d\n",tg1->n,tg2->n,tg3->n,tg4->n);
1180 DAG_traverse(&dt); */
1181 } else { /* no neighbour */
1182 if ((tg1=(triangle *)malloc(sizeof(triangle)))==0) {
1183 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1184 exit(1);
1185 }
1186 if ((tg2=(triangle *)malloc(sizeof(triangle)))==0) {
1187 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1188 exit(1);
1189 }
1190
1191 tg1->p0=p0;tg1->p1=p1;tg1->p2=point_number;tg1->n=ntriangles;
1192 tg2->p0=p0;tg2->p1=point_number;tg2->p2=p3;tg2->n=ntriangles+1;
1193
1194 /* alive, dead triangles */
1195 /* new triangles are alive */
1196 tg1->status=tg2->status=LIVE;
1197 /* old triangle is dead */
1198 tg->status=DEAD;
1199
1200 /* arrange neighbours */
1201 tg1->t0=tg->t2;tg1->t1=tg2->n;tg1->t2=(t0?t0->n:-1);
1202 tg2->t0=tg->t2;tg2->t1=(t3?t3->n:-1);tg2->t2=tg1->n;
1203 /* update old neighbours */
1204 update_neighbours_around(t0,tg->n,tg1->n);
1205 update_neighbours_around(t3,tg->n,tg2->n);
1206
1207 /* insert into RBT */
1208 RBT_insert((void*)tg1, &rt);
1209 RBT_insert((void*)tg2, &rt);
1210 ntriangles+=2;
1211
1212 /* update DAG */
1213 parent_dag_node=tg->dag_p;
1214 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg1);
1215 tg1->dag_p=current_dag_node;
1216 current_dag_node=DAG_insert(&dt, parent_dag_node, (void *)tg2);
1217 tg2->dag_p=current_dag_node;
1218 delaunay_refine(tg1->n,point_number,&rt,&dt);
1219 delaunay_refine(tg2->n,point_number,&rt,&dt);
1220
1221 }
1222 }
1223
1224 }
1225 }
1226
1227 /* free point array */
1228 free(point_array);
1229 /* note: we dont free point_array[i].z values
1230 * because they are stored in the tree */
1231
1232 #ifdef DEBUG
1233 DAG_traverse(&dt);
1234 #endif
1235 /* copy boundary to a temp one */
1236 /*copy_polygon(&bound,&temp_polygon); */
1237 /* create a temporary polygon with all edges */
1238 init_poly(&temp_polygon);
1239 for (i=0; i<nedges; i++){
1240 e.p1 =edge_array[i].p1; e.p2=edge_array[i].p2;
1241 e.type=edge_array[i].type;
1242 #ifdef DEBUG
1243 printf("do_initial_triangulation: building plygon edge (%d,%d)\n",e.p1,e.p2);
1244 #endif
1245 insert_edge_to_poly(&temp_polygon,&e);
1246 }
1247
1248 remove_polygon_mesh_intersections(&temp_polygon);
1249
1250 destroy_polygon(&temp_polygon);
1251
1252 /* remove exterior triangles */
1253 remove_exterior_and_update_triangles(&bound);
1254
1255 /* subdivide triangles */
1256 subdivide_all_triangles();
1257 #ifdef DEBUG
1258 for (i=0; i<M.count;i++) {
1259 printf("triangulate: %d (%lf,%lf) = %lf\n",i,Mx(i,M),My(i,M),Mz(i,M));
1260 }
1261 #endif
1262
1263 }
1264
1265 /* function to generate mesh and solve */
1266 int
solve_problem_from_scratch(const char * filename)1267 solve_problem_from_scratch(const char *filename)
1268 {
1269
1270 /* set the degree of freedom of points */
1271 degree_of_freedom=1;
1272 /* (re) generate mesh */
1273 generate_mesh(filename);
1274
1275 /* solve mesh */
1276 if ( solve_equation!=POISSON && solve_equation!=POISSON_SPARSE ) {
1277 /* number of eigenvectors to be found */
1278 if (requested_degree_of_freedom) {
1279 degree_of_freedom=requested_degree_of_freedom; /* global change */
1280 } else {
1281 degree_of_freedom=3;
1282 }
1283 if((eig_array= (MY_DOUBLE*) realloc((void*)eig_array,(size_t)degree_of_freedom*sizeof(MY_DOUBLE)))==0){
1284 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1285 exit(1);
1286 }
1287 if ( solve_equation==HELMHOLTZ) {
1288 re_number_and_solve_helmholtz();
1289 } else if ((solve_equation == HELMHOLTZ_INHOMO)
1290 || (solve_equation == HELMHOLTZ_FREQ)
1291 || (solve_equation == HELMHOLTZ_BETA)) {
1292 buildup_hash_table_and_equation(M,dt);
1293 }
1294 } else if (solve_equation==POISSON_SPARSE) {
1295 re_number_and_solve_poisson_sparse();
1296 } else { /* POISSON */
1297 re_number_and_solve_poisson();
1298 }
1299
1300 return(0);
1301 }
1302
1303
1304 /* free all allocated variables */
1305 void
free_everything(void)1306 free_everything(void)
1307 {
1308 int i;
1309 /* we dont free anything if M.count==0 */
1310 /* because we assume we have not allocated any memory */
1311 if (M.count) {
1312 /* Mesh */
1313 BIT_free(&M);
1314 /* RBT */
1315 RBT_free(&rt);
1316 /* DAG */
1317 DAG_free(&dt);
1318
1319 /* free edge array */
1320 free(edge_array);
1321 edge_array=0;
1322 /* free boundaries */
1323 for (i=0;i<bound.npoly;i++) {
1324 destroy_polygon(&(bound.poly_array[i]));
1325 }
1326 free(bound.poly_array);
1327 bound.poly_array=0;
1328 bound.npoly=0;
1329 }
1330 /* initialize bitree for 100 points intially */
1331 BIT_init(&M,100);
1332 /* initialize DAG */
1333 DAG_init(&dt,&inclusion,&print_detail);
1334 /* initialize RBT with triangles */
1335 RBT_init(&rt,&equ,<,&print_detail,&delete_triangle);
1336
1337 }
1338
1339 /* function to print help */
1340 static void
print_help(void)1341 print_help(void)
1342 {
1343 printf("usage: "PACKAGE" [options] \n\n");
1344 printf("options are:\n");
1345 printf("-i (FILE) : input file name is 'FILE'\n");
1346 printf("-a (val) : set min area of each triangle to 'val'\n");
1347 printf("-A (val) : set max area of each triangle to 'val'\n");
1348 printf("-h, -? : help, print info\n");
1349 /* printf("-v : be verbose\n");
1350 printf("-b (val) : set max badness of each triangle to 'val'\n");
1351 printf("-s (val) : set splitting priority to 'val'\n");
1352 printf("-l (val) : set max nodes to 'val'\n");
1353 printf("-r (val) : limit iterations to 'val'\n");
1354 printf("-conj : use preconditioned conjugate gradient method to solve\n");
1355 printf("-do : output data file of mesh as \"data.in\"\n"); */
1356 printf("-c (val) : set contour levels to 'val'\n");
1357 printf("-d (contour|grad|fill): display contours,gradient,colourfill by default\n");
1358 printf("-e (poisson|helmholtz|helmin|helmk|helmbeta): type of equation to solve\n");
1359 printf("-s : use sparse solver (works only for Poisson Equation)\n");
1360 printf("-x : do not go into display mode, dump output to file 'FILE.eps'\n");
1361 printf("-V : print version\n");
1362 }
1363
1364 /* function to print version */
1365 static void
print_version(void)1366 print_version(void)
1367 {
1368 printf(PACKAGE" "VERSION", Copyright (C) 2001-2005 Sarod Yatawatta\n");
1369 printf(PACKAGE" comes with ABSOLUTELY NO WARRANTY\n");
1370 printf("This is free software, and you are welcome to redistribute\n");
1371 printf("it under certain conditions given by GNU GENERAL PUBLIC LICENSE\n");
1372 printf("$Id: pdnmesh.c,v 1.57 2005/04/24 05:32:24 sarod Exp $\n");
1373 printf("This is a beta release. Please report bugs to\n");
1374 printf(" <"PACKAGE_BUGREPORT">. Thank you.\n");
1375 }
1376
1377
1378 #ifndef WIN32
1379 void
initialize_global_variables(int argc,char * argv[])1380 initialize_global_variables(int argc, char *argv[])
1381 {
1382
1383 int c;
1384 opterr=0;
1385
1386 contour_levels=40; /* no of contour levels */
1387 current_plotting_contour=0; /* number in z array < degree_of_freedom */
1388 g_badness_limit=2; /* max triangle badness */
1389 g_area_floor=0.001; /* min triangle area */
1390 g_area_ceil=0.005; /* max triangle area */
1391 max_iteration_limit=20;
1392
1393 /* default waveguide parameters */
1394 global_wave_k0=1.0;
1395 global_wave_beta=10.0;
1396
1397 g_gradient_limit=1000;
1398
1399 solve_equation = POISSON; /* POISSON, HELMHOLTZ */
1400
1401 /* default plotting */
1402 plot_mesh=0;
1403 plot_cont=0;
1404 plot_grad=0;
1405 plot_legend=0;
1406 text_mode_only=0;
1407
1408 mesh_already_present=0;
1409
1410 update_status_bar("Starting intialization...");
1411 /* getopt loop */
1412 while ((c=getopt(argc,argv,"i:e:d:a:A:c:xsvhV"))!=-1) {
1413 switch (c) {
1414 case 'i':
1415 if (optarg) {
1416 input_cord_filename=(char*)realloc((void*)input_cord_filename,sizeof(char)*(size_t)(strlen((char*)optarg)+1));
1417 if ( input_cord_filename == 0 ) {
1418 fprintf(stderr,"%s: %d: no free memory",__FILE__,__LINE__);
1419 exit(1);
1420 }
1421 strcpy(input_cord_filename,(char*)optarg);
1422 } else {
1423 input_cord_filename=0;
1424 }
1425 break;
1426 case 'e':
1427 if (!strcmp(optarg,"helmholtz")) {
1428 solve_equation=HELMHOLTZ;
1429 } else if( !strcmp(optarg,"poisson")) {
1430 solve_equation=POISSON;
1431 } else if( !strcmp(optarg,"helmin")) {
1432 solve_equation=HELMHOLTZ_INHOMO;
1433 } else if( !strcmp(optarg,"helmk")) {
1434 solve_equation=HELMHOLTZ_FREQ;
1435 } else if( !strcmp(optarg,"helmbeta")) {
1436 solve_equation=HELMHOLTZ_BETA;
1437 } else {
1438 fprintf(stderr,"Unknown equation type `%s`.\n",optarg);
1439 }
1440 break;
1441 case 'd':
1442 if (!strcmp(optarg,"contour")) {
1443 plot_cont=1;
1444 } else if( !strcmp(optarg,"grad")) {
1445 plot_grad=1;
1446 } else if( !strcmp(optarg,"fill")) {
1447 plot_fill=1;
1448 } else {
1449 fprintf(stderr,"Unknown display option `%s`.\n",optarg);
1450 }
1451 break;
1452 case 'a':
1453 g_area_floor=strtod(optarg,0);
1454 if ( g_area_floor == 0.0f ) {
1455 fprintf(stderr,"Invalid value for area `%s`.\n",optarg);
1456 g_area_floor=0.001;
1457 }
1458 break;
1459 case 'A':
1460 g_area_ceil=strtod(optarg,0);
1461 if ( g_area_ceil== 0.0f ) {
1462 fprintf(stderr,"Invalid value for area `%s`.\n",optarg);
1463 g_area_ceil=0.005;
1464 }
1465 break;
1466 case 'c':
1467 contour_levels=strtol(optarg,0,0);
1468 if ( !contour_levels ) {
1469 fprintf(stderr,"Invalid value for contours `%s`.\n",optarg);
1470 contour_levels=40;
1471 }
1472 break;
1473 case 's':
1474 if(solve_equation==POISSON) solve_equation=POISSON_SPARSE;
1475 break;
1476 case 'x':
1477 text_mode_only=1;
1478 break;
1479 case 'h':
1480 print_help();
1481 exit(0);
1482 case 'V':
1483 print_version();
1484 exit(0);
1485
1486 case '?':
1487 if (isprint(optopt))
1488 fprintf(stderr,"Unknown option `-%c`.\n",optopt);
1489 else
1490 fprintf(stderr,
1491 "Unknown option character `\\x%x`.\n",optopt);
1492 exit(1);
1493 default:
1494 break;
1495 }
1496 }
1497
1498 if (!plot_cont && !plot_grad && !plot_fill) {
1499 plot_cont=1;
1500 }
1501 /* initialize bitree for 100 points intially */
1502 BIT_init(&M,100);
1503 /* initialize DAG */
1504 DAG_init(&dt,&inclusion,&print_detail);
1505 /* initialize RBT with triangles */
1506 RBT_init(&rt,&equ,<,&print_detail,&delete_triangle);
1507
1508
1509
1510 /* set eig_array to null */
1511 if((eig_array= (MY_DOUBLE*) calloc((size_t)degree_of_freedom,sizeof(MY_DOUBLE)))==0){
1512 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1513 exit(1);
1514 }
1515
1516 update_status_bar("Done intialization.");
1517 }
1518 #endif /*WIN32 */
1519 /******************************************/
1520
1521 #ifndef WIN32
1522 int
main(int argc,char ** argv)1523 main(int argc, char **argv)
1524 {
1525 /* gtkglext */
1526 GtkWidget *window;
1527 GdkGLConfig *glconfig;
1528 char *output_filename;
1529
1530 /* use getopt to parse commandline */
1531 initialize_global_variables(argc,argv);
1532
1533 if ( input_cord_filename ) {
1534 solve_problem_from_scratch(input_cord_filename);
1535 }
1536 /*free_everything(); */
1537
1538 if(text_mode_only && input_cord_filename) {
1539 /* allocate memory for output filename */
1540 if((output_filename= (char*) malloc((strlen(input_cord_filename)+5)*sizeof(char)))==0){
1541 fprintf(stderr, "%s: %d: no free memory\n", __FILE__,__LINE__);
1542 exit(1);
1543 }
1544 strcpy(output_filename,input_cord_filename);
1545 strcat(output_filename,".eps");
1546 /* save output file */
1547 if (plot_cont) {
1548 print_contour_all(output_filename,1,M);
1549 }
1550 if (plot_grad) {
1551 print_gradient(output_filename,M);
1552 }
1553
1554 free(output_filename);
1555 return(0);
1556 }
1557 /* intialiaze gtk */
1558 gtk_init(&argc,&argv);
1559 /* intialize GtkGlExt */
1560 gtk_gl_init(&argc,&argv);
1561 /* configure OpenGL framebuffer */
1562 glconfig=configure_gl();
1563
1564 /* create and show the application window */
1565 window=create_window(glconfig);
1566 gtk_widget_show(window);
1567 gtk_main();
1568 return(0);
1569 }
1570 #endif/*WIN32 */
1571