1 /* Copyright (C) 1992-1998 The Geometry Center
2 * Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips
3 *
4 * This file is part of Geomview.
5 *
6 * Geomview is free software; you can redistribute it and/or modify it
7 * under the terms of the GNU Lesser General Public License as published
8 * by the Free Software Foundation; either version 2, or (at your option)
9 * any later version.
10 *
11 * Geomview is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * Lesser General Public License for more details.
15 *
16 * You should have received a copy of the GNU Lesser General Public
17 * License along with Geomview; see the file COPYING. If not, write
18 * to the Free Software Foundation, 675 Mass Ave, Cambridge, MA 02139,
19 * USA, or visit http://www.gnu.org.
20 */
21
22 #if HAVE_CONFIG_H
23 # include "config.h"
24 #endif
25
26 #if 0
27 static char copyright[] = "Copyright (C) 1992-1998 The Geometry Center\n\
28 Copyright (C) 1998-2000 Stuart Levy, Tamara Munzner, Mark Phillips";
29 #endif
30
31 /* Comments on numerical issues: Dirichlet() uses two constants,
32 namely VERTEX_EPSILON and MATRIX_EPSILON, to decide when two
33 real numbers are equal. Setting the values of these constants
34 is tricky. In situations where the polyhedron has very tiny
35 faces and the initial list of generators is sufficiently
36 precise (e.g. (40, 1) surgery on the figure eight knot) the
37 constants should be set fairly small. In situations where the
38 polyhedron has no small faces and the generators are very
39 imprecise (e.g. (1,-1)(-5,1)(-5,1) surgery on the Borromean
40 rings, which gives the manifold of volume 0.94...) the
41 constants must be fairly large. (The more complicated
42 Dehn fillings give less precise generators not because the
43 original solution to the gluing equations is less accurate,
44 but because the original set of generators is far from the
45 canonical set, and a lot of roundoff error accumulates in
46 the basepoint finding routine.) */
47
48 /* #define PRECISE_GENERATORS 1 for high precision generators
49 and tiny faces, or #define PRECISE_GENERATORS 0 for low
50 precision generators and large faces. */
51 #define PRECISE_GENERATORS 0
52
53 /* Dirichlet() uses an algorithm based on checking that
54 that corresponding faces do in fact match up. Specifically,
55 for each face it checks that (the image of) its matching
56 face does not extend beyond the perimeter of the given face.
57 It checks that the matching face does not
58 extend beyond any Dirichlet plane determined by a group
59 element (alpha)(beta), where alpha is the group element
60 associated with the matching face, and beta is the group
61 element associated with one of the faces bordering the given
62 face. When the faces do match up
63 correctly the program can verify the fact very quickly, because
64 in the typical case where the edges are of order three the
65 group element (alpha)(beta) is exactly the group element
66 associated with a face bordering on the matching
67 face. So the program can traverse the original face
68 counterclockwise while simultaneously traversing the matching
69 face clockwise and noting that group elements correspond
70 correctly. If a group element fails to correspond, then the
71 program must check that no vertices of the matching face
72 lie beyond the Dirichlet plane corresponding to (alpha)(beta).
73 But even this is pretty fast--it just needs to evaluate some
74 inner products and check that none are positive. In the event
75 that some vertices do lie beyond the Dirichlet plane
76 corresponding to (alpha)(beta), the program knows to add this
77 plane (and its inverse) to the collection of faces of the
78 polyhedron. Thus the program adds only faces that are
79 actually needed. This allows it to quickly find Dirichlet
80 domains whose vertices (and face planes) are too far from the
81 origin to be found by enumerating all group elements
82 at that distance (e.g. (40,1) Dehn filling on a knot
83 complement).
84
85 The results of this algorithm can be made rigorous by checking
86 that the sum of the dihedral angles around each edge add up
87 to 2pi and that the face pairings are sufficient to generate
88 the original group, but I haven't written the code to do this
89 yet. (In fact, the condition that the face pairings generate
90 the original group might be satisfied automatically given
91 that we start with a generating set that includes all
92 group elements up to some norm. Perhaps there's also a
93 way to prove that the edge angles must add up correctly.)
94 */
95
96 /* Note: Dirichlet() does not require inverses of generators */
97 /* to be present. It provides any inverses which may be */
98 /* missing. */
99
100 #include "options.h"
101 #include "projective.h"
102 #include "complex.h"
103 #include "winged_edge.h"
104 #include "3d.h"
105 #include "dgflag.h"
106 #include <stdio.h>
107 #include "extern.h"
108 #include "discgrpP.h"
109
110 #define free32(ptr) free(ptr)
111 #define malloc32(size) malloc(size)
112 #if PRECISE_GENERATORS
113 #define VERTEX_EPSILON (1e9 * HARDWARE_PRECISION)
114 #define MATRIX_EPSILON (1e7 * HARDWARE_PRECISION)
115 #else
116 #define VERTEX_EPSILON 1e-3
117 #define MATRIX_EPSILON 1e-5
118 #endif
119
120 #define MIN_DL_FILL_TONE 0
121 #define MAX_DL_FILL_TONE 64
122
123 typedef double vector[4];
124
125
126 /*static void convert_generators(sl2c_matrix *m, proj_matrix **mm, int foo);*/
127 static void make_cube(WEpolyhedron *poly);
128 /*static void make_hypercube(WEpolyhedron *poly);*/
129 static void initialize_polyhedron(WEpolyhedron *poly, proj_matrix *m, int n);
130 static int find_Dirichlet_domain(WEpolyhedron *poly);
131 static int check_face(WEpolyhedron *poly, WEface *face);
132 static int all_dirty_faces_unmatched(WEpolyhedron *poly);
133 static int unsophisticated_search(WEpolyhedron *poly);
134 static int add_element(WEpolyhedron *poly, proj_matrix m);
135 static void slice_off_cusps(WEpolyhedron *poly);
136 static int add_face(WEpolyhedron *poly, proj_matrix m, WEface *face);
137 static void cut_edges(WEpolyhedron *poly);
138 static void adjust_f_e_ptrs(WEpolyhedron *poly);
139 static void cut_faces(WEpolyhedron *poly, WEface *face);
140 static void remove_dead_edges(WEpolyhedron *poly);
141 static void remove_dead_vertices(WEpolyhedron *poly);
142
143 static void number_faces(WEpolyhedron *poly);
144 /*static void read_vertices(WEface *face, double (*foo)[3]);*/
145 /*static void print_statistics(WEpolyhedron *poly);*/
146 static void print_poly(WEpolyhedron *poly);
147 static void print_vef(WEpolyhedron *poly);
148 /*static void print_vertex_distances(WEpolyhedron *poly);*/
149 static void print_vertices(WEpolyhedron *poly);
150 /*static void print_edge_lengths(WEpolyhedron *poly);*/
151 /*static void print_face_distances(WEpolyhedron *poly);*/
152 /*static void saveOOGL(WEpolyhedron *poly);*/
153 /*static void free_polyhedron(WEpolyhedron *poly);*/
154 static void roundoff_message(char *message);
155
156 static int matrix_epsilon_message_given,
157 vertex_epsilon_message_given;
158
159 static point origin;
160 static int metric, debug;
161 static WEpolyhedron the_polyhedron;
162
163 /* we'd use the geometry library routines here if they were double precision */
164 extern double DHPt3Dot3();
165 extern double DHPt3Dot();
166
167 static WEedge edgeheap[1024];
168 static int heapcount = 0;
169
170 int proj_same_matrix();
171
172 void
do_weeks_code(WEpolyhedron ** wepolyhedron,point oo,proj_matrix * gen_list,int n,int met,int slice)173 do_weeks_code(WEpolyhedron **wepolyhedron, point oo, proj_matrix *gen_list,
174 int n, int met, int slice)
175 {
176 int i;
177 for (i=0; i<4; ++i) origin[i] = oo[i];
178
179 heapcount = 0;
180 metric = met;
181 *wepolyhedron = &the_polyhedron;
182 /*
183 if (metric == DG_SPHERICAL) make_hypercube(*wepolyhedron);
184 else make_cube(*wepolyhedron);
185 */
186 make_cube(*wepolyhedron);
187 initialize_polyhedron(*wepolyhedron, gen_list, n);
188 if (find_Dirichlet_domain(*wepolyhedron) == 0) *wepolyhedron = NULL;
189 else {
190 if (debug == 2) print_poly(*wepolyhedron);
191 number_faces(*wepolyhedron);
192 }
193 if (slice && (metric == DG_HYPERBOLIC)) slice_off_cusps(*wepolyhedron);
194 }
195
196 #define MAGIC_SCALE .99
197 static void
slice_off_cusps(poly)198 slice_off_cusps(poly)
199 WEpolyhedron *poly;
200 {
201 WEvertex *vptr, *vlist;
202 Transform tlate;
203 proj_matrix dtlate;
204 HPoint3 tlatept;
205 int i,j;
206
207 /* make a copy of the vertex list as it stands now */
208 vlist = (WEvertex *) malloc32(poly->num_vertices * sizeof(WEvertex));
209 vptr = poly->vertex_list;
210 i = 0;
211 do {
212 vlist[i] = *vptr;
213 vlist[i].next = &vlist[i+1];
214 vptr = vptr->next;
215 i++;
216 } while (vptr != NULL);
217 vlist[i-1].next = NULL;
218
219 /* now cycle through these vertices */
220 vptr = vlist;
221 do {
222 vptr->ideal = ( (DHPt3Dot(vptr->x, vptr->x, metric)) >= -.00005) ? 1 : 0;
223 if (vptr->ideal) {
224 tlatept.x = vptr->x[0] * MAGIC_SCALE;
225 tlatept.y = vptr->x[1] * MAGIC_SCALE;
226 tlatept.z = vptr->x[2] * MAGIC_SCALE;
227 tlatept.w = vptr->x[3];
228 Tm3HypTranslateOrigin(tlate, &tlatept);
229 for (i=0;i<4;++i)
230 for (j=0;j<4;++j)
231 dtlate[j][i] = tlate[i][j];
232 add_element(poly, dtlate);
233 }
234 vptr = vptr->next;
235 } while (vptr != NULL);
236
237 }
238
239 #if 0
240 static void convert_generators(sl2c_generators, proj_generators_ptr, num_generators)
241 sl2c_matrix *sl2c_generators;
242 proj_matrix **proj_generators_ptr;
243 int num_generators;
244 {
245 int i;
246
247 *proj_generators_ptr = (proj_matrix *) malloc32((MyInt32) num_generators * sizeof(proj_matrix));
248
249 for (i=num_generators; --i>=0; )
250 sl2c_to_proj(sl2c_generators[i], (*proj_generators_ptr)[i]);
251
252 return;
253 }
254 #endif
255
256 #if 0
257 static void make_hypercube(polyhedron)
258 WEpolyhedron *polyhedron;
259 {
260 int i;
261 WEvertex *initial_vertices[16];
262 WEedge *initial_edges[24];
263 WEface *initial_faces[12];
264 static int edata[24][8] = {
265 {0, 4, 8, 4, 9, 6, 1,0},
266 {2, 6, 4, 16, 6, 11, 0,8},
267 {1, 5, 20, 8, 7, 9, 5,1},
268 {3, 11, 10, 5, 22, 17, 7,6},
269 {0, 2, 0, 8, 1, 10, 0,2},
270 {1, 3, 8, 20, 10, 3, 2,6},
271 {4, 6, 12, 0, 11, 1, 3,0},
272 {5, 13, 2, 9, 14, 21, 5, 4},
273 {0, 1, 4, 0, 5, 2, 2,1},
274 {4, 5, 0, 12, 2, 7, 1,4},
275 {2, 3, 16, 4, 3, 5, 7,2},
276 {6, 14, 6, 1, 18, 13, 3,8},
277
278 {4, 4+8, 9, 6, 9+12, 6+12, 4,3},
279 {2+8, 6+8, 4+12, 10+12, 11, 11+12, 8,10},
280 {1+8, 5+8, 5+12, 8+12, 7+12, 7, 9,5},
281 {3+8, 7+8, 10+12, 5+12, 11+12, 7+12, 10,9},
282 {2, 2+8, 1, 10, 1+12, 10+12, 8,7},
283 {1+8, 3+8, 8+12, 2+12, 3, 3+12, 6,9},
284 {4+8, 6+8, 9+12, 0+12, 11+12, 11, 11,3},
285 {5+8, 7+8, 2+12, 9+12, 3+12, 11+12, 9,11},
286 {1, 1+8, 5, 2, 5+12, 2+12, 6,5},
287 {4+8, 5+8, 0+12, 6+12, 7, 7+12, 4,11},
288 {2+8, 3+8, 1+12, 4+12, 3+12, 3, 10,7},
289 {6+8, 7+8, 6+12, 1+12, 7+12, 3+12, 11,10}};
290 static int fdata[12] = {0, 8,4,6,9,2,5,10,1,15,23,19};
291
292
293 polyhedron->num_vertices = 16;
294 polyhedron->num_edges = 24;
295 polyhedron->num_faces = 12;
296
297 /*
298 vt = (WEvertex *) malloc32( 16 * sizeof(WEvertex));
299 et = (WEedge *) malloc32( 24 * sizeof(WEedge));
300 ft = (WEface *) malloc32( 12 * sizeof(WEface));
301 */
302
303 for (i=16; --i>=0; )
304 initial_vertices[i] = (WEvertex *) malloc32((MyInt32) sizeof(WEvertex));
305 for (i=24; --i>=0; )
306 initial_edges[i] = (WEedge *) malloc32((MyInt32) sizeof(WEedge));
307 for (i=12; --i>=0; )
308 initial_faces[i] = (WEface *) malloc32((MyInt32) sizeof(WEface));
309
310 /*
311 for (i=0; i<16; ++i)
312 initial_vertices[i] = &vt[i];
313 for (i=0; i<24; ++i )
314 initial_edges[i] = &et[i];
315 for (i=0; i <12; ++i )
316 initial_faces[i] = &ft[i];
317 */
318
319 polyhedron->vertex_list = initial_vertices[0];
320 polyhedron->edge_list = initial_edges[0];
321 polyhedron->face_list = initial_faces[0];
322
323 polyhedron->dirty0.nxt = initial_faces[0];
324 polyhedron->dirty0.prv = NULL; /* should be unnecessary */
325 polyhedron->dirty1.nxt = NULL; /* should be unnecessary */
326 polyhedron->dirty1.prv = initial_faces[11];
327
328 polyhedron->clean0.nxt = &polyhedron->clean1;
329 polyhedron->clean0.prv = NULL; /* should be unnecessary */
330 polyhedron->clean1.nxt = NULL; /* should be unnecessary */
331 polyhedron->clean1.prv = &polyhedron->clean0;
332
333 for (i=0; i<16; ++i ) {
334 initial_vertices[i]->x[0] = (i & 8) ? 1.0 : -1.0;
335 initial_vertices[i]->x[1] = (i & 4) ? 1.0 : -1.0;
336 initial_vertices[i]->x[2] = (i & 2) ? 1.0 : -1.0;
337 initial_vertices[i]->x[3] = (i & 1) ? 1.0 : -1.0;
338 initial_vertices[i]->next = initial_vertices[i+1];
339 }
340 initial_vertices[15]->next = NULL;
341 /* delete vertices # 7 and 8 */
342 initial_vertices[6]->next = initial_vertices[9];
343 polyhedron->num_vertices = 14;
344
345 for (i=0; i< 24; ++i ) {
346 initial_edges[i]->v0 = initial_vertices[edata[i][0]];
347 initial_edges[i]->v1 = initial_vertices[edata[i][1]];
348 initial_edges[i]->e0L = initial_edges[edata[i][2]];
349 initial_edges[i]->e0R = initial_edges[edata[i][3]];
350 initial_edges[i]->e1L = initial_edges[edata[i][4]];
351 initial_edges[i]->e1R = initial_edges[edata[i][5]];
352 initial_edges[i]->fL = initial_faces[edata[i][6]];
353 initial_edges[i]->fR = initial_faces[edata[i][7]];
354 initial_edges[i]->next = initial_edges[i+1];
355 }
356 initial_edges[23]->next = NULL;
357
358 for (i=0; i<12; ++i) {
359 initial_faces[i]->order = 4;
360 initial_faces[i]->fill_tone = -2;
361 initial_faces[i]->some_edge = initial_edges[fdata[i]];
362 initial_faces[i]->inverse = NULL;
363 initial_faces[i]->next = initial_faces[i+1];
364 initial_faces[i]->prv = initial_faces[i-1];
365 initial_faces[i]->nxt = initial_faces[i+1];
366 }
367 initial_faces[11]->next = NULL;
368 initial_faces[0]->prv = &polyhedron->dirty0;
369 initial_faces[11]->nxt = &polyhedron->dirty1;
370
371 return;
372 }
373 #endif
374
make_cube(polyhedron)375 static void make_cube(polyhedron)
376 WEpolyhedron *polyhedron;
377 {
378 int i;
379 WEvertex *initial_vertices[8];
380 WEedge *initial_edges[12];
381 WEface *initial_faces[6];
382 static int edata[12][8] = {
383 {0, 4, 8, 4, 9, 6, 2, 4},
384 {2, 6, 4, 10, 6, 11, 4, 3},
385 {1, 5, 5, 8, 7, 9, 5, 2},
386 {3, 7, 10, 5, 11, 7, 3, 5},
387 {0, 2, 0, 8, 1, 10, 4, 0},
388 {1, 3, 8, 2, 10, 3, 0, 5},
389 {4, 6, 9, 0, 11, 1, 1, 4},
390 {5, 7, 2, 9, 3, 11, 5, 1},
391 {0, 1, 4, 0, 5, 2, 0, 2},
392 {4, 5, 0, 6, 2, 7, 2, 1},
393 {2, 3, 1, 4, 3, 5, 3, 0},
394 {6, 7, 6, 1, 7, 3, 1, 3}};
395 static int fdata[6] = {4, 6, 0, 1, 0, 2};
396
397
398 polyhedron->num_vertices = 8;
399 polyhedron->num_edges = 12;
400 polyhedron->num_faces = 6;
401
402 for (i=8; --i>=0; )
403 initial_vertices[i] = (WEvertex *) malloc32((MyInt32) sizeof(WEvertex));
404 for (i=12; --i>=0; )
405 initial_edges[i] = (WEedge *) malloc32((MyInt32) sizeof(WEedge));
406 for (i=6; --i>=0; )
407 initial_faces[i] = (WEface *) malloc32((MyInt32) sizeof(WEface));
408
409 polyhedron->vertex_list = initial_vertices[0];
410 polyhedron->edge_list = initial_edges[0];
411 polyhedron->face_list = initial_faces[0];
412
413 polyhedron->dirty0.nxt = initial_faces[0];
414 polyhedron->dirty0.prv = NULL; /* should be unnecessary */
415 polyhedron->dirty1.nxt = NULL; /* should be unnecessary */
416 polyhedron->dirty1.prv = initial_faces[5];
417
418 polyhedron->clean0.nxt = &polyhedron->clean1;
419 polyhedron->clean0.prv = NULL; /* should be unnecessary */
420 polyhedron->clean1.nxt = NULL; /* should be unnecessary */
421 polyhedron->clean1.prv = &polyhedron->clean0;
422
423 for (i=8; --i>=0; ) {
424 initial_vertices[i]->x[0] = (i & 4) ? 17.0 : -17.0;
425 initial_vertices[i]->x[1] = (i & 2) ? 17.0 : -17.0;
426 initial_vertices[i]->x[2] = (i & 1) ? 17.0 : -17.0;
427 /*
428 if (metric & DG_SPHERICAL)
429 initial_vertices[i]->x[3] = 0.0;
430 else
431 */
432 initial_vertices[i]->x[3] = 1.0;
433 initial_vertices[i]->next = initial_vertices[i+1];
434 }
435 initial_vertices[7]->next = NULL;
436
437 for (i=12; --i>=0; ) {
438 initial_edges[i]->v0 = initial_vertices[edata[i][0]];
439 initial_edges[i]->v1 = initial_vertices[edata[i][1]];
440 initial_edges[i]->e0L = initial_edges[edata[i][2]];
441 initial_edges[i]->e0R = initial_edges[edata[i][3]];
442 initial_edges[i]->e1L = initial_edges[edata[i][4]];
443 initial_edges[i]->e1R = initial_edges[edata[i][5]];
444 initial_edges[i]->fL = initial_faces[edata[i][6]];
445 initial_edges[i]->fR = initial_faces[edata[i][7]];
446 initial_edges[i]->next = initial_edges[i+1];
447 }
448 initial_edges[11]->next = NULL;
449
450 for (i=6; --i>=0; ) {
451 initial_faces[i]->order = 4;
452 initial_faces[i]->fill_tone = -2;
453 initial_faces[i]->some_edge = initial_edges[fdata[i]];
454 initial_faces[i]->inverse = NULL;
455 if (i < 5) {
456 initial_faces[i]->next = initial_faces[i+1];
457 }
458 if (i > 0) {
459 initial_faces[i]->prv = initial_faces[i-1];
460 }
461 initial_faces[i]->nxt = initial_faces[i+1];
462 }
463 initial_faces[5]->next = NULL;
464 initial_faces[0]->prv = &polyhedron->dirty0;
465 initial_faces[5]->nxt = &polyhedron->dirty1;
466
467 return;
468 }
469
470
initialize_polyhedron(polyhedron,proj_generators,num_generators)471 static void initialize_polyhedron(polyhedron, proj_generators, num_generators)
472 WEpolyhedron *polyhedron;
473 proj_matrix *proj_generators;
474 int num_generators;
475 {
476 int i;
477 WEface *face;
478
479 for (i=num_generators; --i>=0; )
480 {
481 add_element(polyhedron, proj_generators[i]);
482 if (debug == 2) print_poly(polyhedron);
483 }
484
485 /* make sure no faces of the original cube remain */
486 for (face=polyhedron->face_list; face; face=face->next)
487 if (debug) if (face->fill_tone == -2) {
488 fprintf(stderr, "A face of the original cube is inside the polyhedron\n");
489 fprintf(stderr, "determined by the initial generators. This program\n");
490 fprintf(stderr, "could be modified to deal with this situation, but\n");
491 fprintf(stderr, "it's not ready yet.\n");
492 /* the solution would be to tile to some small */
493 /* radius (or increase the size of the cube!) */
494 return;
495 }
496
497 return;
498 }
499
check_face(polyhedron,face)500 static int check_face(polyhedron, face)
501 WEpolyhedron *polyhedron;
502 WEface *face;
503 {
504 /* We want to see whether the image of the matching face */
505 /* is entirely contained within the given face. To do this */
506 /* we'll examine the edges of the given face one at a time */
507 /* and make sure the image of the matching face doesn't */
508 /* extend past any edge. We do this as follows. Let alpha */
509 /* be the inverse of the group element associated with the */
510 /* given face, and beta be the group element associated */
511 /* with one of its neighboring faces. We want to make sure */
512 /* that none of the vertices of the matching face extend */
513 /* beyond the Dirichlet plane determined by (alpha)(beta). */
514 /* This will be true for all beta iff the image of the */
515 /* matching face is contained within the given face. */
516 /* */
517 /* We want this routine to run correctly in all cases and */
518 /* to run quickly in the typical case. In the typical case */
519 /* the faces will already match exactly, and the edges will */
520 /* all be of order three, so (alpha)(beta) will in fact */
521 /* be the group element associated with one of the faces */
522 /* bordering on the matching face. So in this case we */
523 /* simply traverse the given face in the counterclockwise */
524 /* while simultaneously traversing the matching face in the */
525 /* clockwise direction, and check that each (alpha)(beta) */
526 /* coincides with a gamma bordering the matching face. If */
527 /* this fails (as will occur when the faces do not coincide */
528 /* or the edges are not of order 3) then our fallback plan */
529 /* is to check that all vertices of the matching face lie */
530 /* on the correct side of the Dirichlet plane determined by */
531 /* (alpha)(beta). */
532 /* */
533 /* (Note: there may be a better way to structure the */
534 /* algorithm so that we check edges instead of faces. For */
535 /* a typical edge (of order three) we want to check that */
536 /* the product of the three surrounding group elements is */
537 /* the identity. But that will wait for another day.) */
538
539 WEface *match_face;
540 WEedge *edge,
541 *edge0,
542 *placeholder;
543 WEvertex *vertex;
544 double (*alpha)[4], /* these are proj_matrices with */
545 (*beta)[4], /* no storage attached */
546 (*gamma)[4];
547 vector planne;
548 proj_matrix alphabeta;
549 int alphabeta_found;
550 point gorigin;
551
552 /* If the given face doesn't have an active matching face */
553 /* (i.e. its mate lies entirely outside the polyhedron) */
554 /* then we put it at the end of the dirty list and hope it */
555 /* goes away. */
556 /* Dec. 8, 92: don't bother with putting the original faces of */
557 /* the cube on the dirty list. */
558 if (face->inverse == NULL) {
559 if (face->fill_tone != -2) {
560 face->nxt = &polyhedron->dirty1;
561 face->prv = polyhedron->dirty1.prv;
562 polyhedron->dirty1.prv->nxt = face;
563 polyhedron->dirty1.prv = face;
564 polyhedron->pending0.nxt = &polyhedron->pending1;
565 polyhedron->pending1.prv = &polyhedron->pending0;
566 }
567
568 /* Be sure the dirty queue contains at least one */
569 /* face with an active inverse face. If it doesn't */
570 /* call a routine to try to create some new faces. */
571 /* If that fails, print a message and exit. */
572 while (all_dirty_faces_unmatched(polyhedron))
573 {
574 if (debug)
575 fprintf(stderr, "searching for more group elements\n");
576 if ( ! unsophisticated_search(polyhedron)) {
577 if (debug) fprintf(stderr, "The dirty list in Dirichlet.c contains only unmatched faces.\n");
578 return(0);
579 }
580 }
581 return(1);
582 }
583
584 match_face = face->inverse;
585 alpha = match_face->group_element;
586
587 /* The value of edge0 is maintained from one use of the */
588 /* inner loop to the next to improve efficiency in the */
589 /* case where the faces match exactly. */
590 edge0 = match_face->some_edge;
591
592 /* traverse the face */
593 edge = face->some_edge;
594 do {
595 if (edge->fL == face) { /* edge points counterclockwise */
596 beta = edge->fR->group_element;
597 edge = edge->e1L;
598 }
599 else { /* edge points clockwise */
600 beta = edge->fL->group_element;
601 edge = edge->e0R;
602 }
603
604 proj_mult(alpha, beta, alphabeta);
605
606 /* is alphabeta the matrix associated with a */
607 /* neighbor of match_face? */
608 alphabeta_found = 0;
609 placeholder = edge0;
610 do {
611 /* edge points counterclockwise */
612 if (edge0->fL == match_face) {
613 gamma = edge0->fR->group_element;
614 edge0 = edge0->e0L; /* traverse CLOCKWISE */
615 }
616 else { /* edge points clockwise */
617 gamma = edge0->fL->group_element;
618 edge0 = edge0->e1R; /* traverse CLOCKWISE */
619 }
620
621 if (proj_same_matrix(alphabeta, gamma)) {
622 alphabeta_found = 1;
623 break;
624 }
625 }
626 while (edge0 != placeholder);
627
628 /* If alphabeta_found == 1 we know the image of the */
629 /* match_face doesn't extend beyond this edge. If */
630 /* alphabeta_found == 0 we need to investigate */
631 /* further by checking whether the vertices of the */
632 /* match_face all lie on the correct side of the */
633 /* Dirichlet plane determined by alphabeta. */
634 if ( ! alphabeta_found) {
635 /* see comments in add_face() */
636 /* first get the image of origin */
637 matvecmul4(alphabeta, origin, gorigin);
638 DHPt3PerpBisect( origin, gorigin, planne,metric);
639
640 do {
641 if (edge0->fL == match_face) { /* CCL */
642 vertex = edge0->v0;
643 edge0 = edge0->e1L;
644 }
645 else { /* CL */
646 vertex = edge0->v1;
647 edge0 = edge0->e0R;
648 }
649
650 if (DHPt3Dot(planne, vertex->x,metric) > VERTEX_EPSILON) {
651 /* add the new group element */
652 add_element(polyhedron, alphabeta);
653
654 /* if face has been modified, go home */
655 if (polyhedron->pending0.nxt != face)
656 return(1);
657
658 /* otherwise if match_face still exists, */
659 /* keep checking face */
660 if (face->inverse) {
661 edge0 = match_face->some_edge;
662 break;
663 }
664
665 /* otherwise put face on the dirty list and */
666 /* go home */
667 face->nxt = &polyhedron->dirty1;
668 face->prv = polyhedron->dirty1.prv;
669 polyhedron->dirty1.prv->nxt = face;
670 polyhedron->dirty1.prv = face;
671 polyhedron->pending0.nxt = &polyhedron->pending1;
672 polyhedron->pending1.prv = &polyhedron->pending0;
673 return(1);
674 }
675 }
676 while (edge0 != placeholder);
677 }
678 }
679 while (edge != face->some_edge);
680
681 return(1);
682 }
683
find_Dirichlet_domain(polyhedron)684 static int find_Dirichlet_domain(polyhedron)
685 WEpolyhedron *polyhedron;
686 {
687 /* By the time this routine is called the faces */
688 /* corresponding to the initial generators should all */
689 /* be on the dirty list, and the clean list should be */
690 /* empty. The faces are examined one at a time to see */
691 /* whether (the image of) the matching face is entirely */
692 /* contained within the given face. If it is, the face */
693 /* is put on the clean list. If it's not, a new face */
694 /* is added to the polyhedron. New faces, and old */
695 /* faces which have been cut, are put on the dirty list.*/
696
697 WEface *face;
698
699 while (polyhedron->dirty0.nxt != &polyhedron->dirty1) {
700 /* pull a face off the dirty list */
701 face = polyhedron->dirty0.nxt;
702 polyhedron->dirty0.nxt = face->nxt;
703 face->nxt->prv = &polyhedron->dirty0;
704
705 /* put the face on the pending list so that if */
706 /* add_face() feels the need to put it on the */
707 /* dirty list, it will have a list to remove */
708 /* it from */
709 face->nxt = &polyhedron->pending1;
710 face->prv = &polyhedron->pending0;
711 polyhedron->pending0.nxt = face;
712 polyhedron->pending1.prv = face;
713
714 /* check the face */
715 if (check_face(polyhedron, face) == 0)
716 /* don't worry if original faces hang around */
717 if (face->fill_tone != -2) return 0;
718
719 /* if the face is still on the pending list, */
720 /* move it to the clean list */
721 if (polyhedron->pending0.nxt == face) {
722 face->nxt = polyhedron->clean0.nxt;
723 face->prv = &polyhedron->clean0;
724 polyhedron->clean0.nxt->prv = face;
725 polyhedron->clean0.nxt = face;
726 }
727 }
728
729 return 1;
730 }
731
732
733
734
all_dirty_faces_unmatched(polyhedron)735 static int all_dirty_faces_unmatched(polyhedron)
736 WEpolyhedron *polyhedron;
737 {
738 WEface *face;
739 int realdirty = 0;
740
741 /* If at least one face on the dirty list has a */
742 /* matching face, or the dirty list is empty, */
743 /* return 0. Otherwise return 1. */
744 /* Make a change to this routine to ignore faces of the original cube */
745
746 if (polyhedron->dirty0.next == &polyhedron->dirty1)
747 return(0); /* shouldn't occur in practice */
748
749 for (face = polyhedron->dirty0.nxt; face != &polyhedron->dirty1; face = face->nxt)
750 {
751 if (face->fill_tone != -2) /* not face of original cube */
752 {
753 if ( face->inverse) return(0);
754 else realdirty = 1; /* there exist unmatched, real faces */
755 }
756 }
757
758 return(realdirty);
759 }
760
761
762 /* unsophisticated_search() returns 1 if it */
763 /* manages to add a new face, 0 if it doesn't. */
764
unsophisticated_search(polyhedron)765 static int unsophisticated_search(polyhedron)
766 WEpolyhedron *polyhedron;
767 {
768 /* This routine serves as a backup for the more efficient */
769 /* algorithm which is normally used. On rare occasions */
770 /* (e.g. for wh_left) the usual algorithm ends up with only */
771 /* unmatched faces on the dirty list, so it cannot proceed. */
772 /* This routine tries looking at all products of two */
773 /* elements in an effort to scare up a new face. */
774
775 WEface *face0,
776 *face1;
777 proj_matrix alpha;
778
779 for (face0=polyhedron->face_list; face0; face0=face0->next)
780 for (face1=polyhedron->face_list; face1; face1=face1->next) {
781 proj_mult(face0->group_element, face1->group_element, alpha);
782 if (add_element(polyhedron, alpha))
783 return(1);
784 }
785 return(0);
786 }
787
788
789
790
791 /* Faces are created in matched pairs (corresponding to a */
792 /* group element and its inverse) by add_element(). If one */
793 /* element of a pair is killed, the inverse pointer of its */
794 /* mate is set to NULL. */
795
796 /* add_element() returns 1 if it adds a face, 0 if it doesn't */
797
add_element(polyhedron,m0)798 static int add_element(polyhedron, m0)
799 WEpolyhedron *polyhedron;
800 proj_matrix m0;
801 {
802 proj_matrix m1;
803 WEface *new_face0,
804 *new_face1;
805 int result0,
806 result1,
807 order2 = 0;
808
809 /* compute the inverse matrix */
810 proj_invert(m0, m1);
811 if (proj_same_matrix(m0, m1)) order2 = 1;
812
813 /* create the new faces */
814 new_face0 = (WEface *) malloc32((MyInt32) sizeof(WEface));
815 new_face1 = (WEface *) malloc32((MyInt32) sizeof(WEface));
816
817 /* set their inverse pointers */
818 new_face0->inverse = new_face1;
819 new_face1->inverse = new_face0;
820
821 /* attempt to add the faces */
822 /* If for any reason one (or both) is unnecessary, */
823 /* add_face() will free the one and adjust the */
824 /* inverse pointer of the other. */
825 if (order2)
826 {
827 new_face0->inverse = new_face0;
828 result0 = add_face(polyhedron, m0, new_face0);
829 return(result0 );
830 }
831 else
832 {
833 result0 = add_face(polyhedron, m0, new_face0);
834 result1 = add_face(polyhedron, m1, new_face1);
835 return(result0 || result1);
836 }
837
838 }
839
840
841 /* add_face() returns 1 if it adds a face, 0 if it doesn't */
842
add_face(polyhedron,matrix,new_face)843 static int add_face(polyhedron, matrix, new_face)
844 WEpolyhedron *polyhedron;
845 proj_matrix matrix;
846 WEface *new_face;
847 {
848 vector planne, gorigin;
849 WEvertex *vertex;
850 int face_is_needed;
851
852 /* get the normal (in the Minkowski metric) */
853 /* to the Dirichlet hyperplane */
854
855 /* first get the image of origin */
856 matvecmul4(matrix,origin, gorigin);
857 DHPt3PerpBisect(origin, gorigin, planne,metric);
858 /* Compute the "distance" from each vertex to the */
859 /* Dirichlet plane. Count the number of vertices */
860 /* on or beyond the plane. Vertices within */
861 /* VERTEX_EPSILON of the plane are assumed to lie */
862 /* on the plane. */
863 face_is_needed = 0;
864 for (vertex=polyhedron->vertex_list; vertex; vertex=vertex->next) {
865 vertex->dist = DHPt3Dot( planne, vertex->x, metric);
866 if (vertex->dist > VERTEX_EPSILON)
867 face_is_needed = 1;
868 else if (vertex->dist > - VERTEX_EPSILON) {
869 if (fabs(vertex->dist) > 1e-2 * VERTEX_EPSILON && ! vertex_epsilon_message_given) {
870 if (debug) roundoff_message("VERTEX_EPSILON");
871 vertex_epsilon_message_given = 1;
872 }
873 vertex->dist = 0.0;
874 }
875 }
876
877 if ( ! face_is_needed) {
878 if (new_face->inverse)
879 new_face->inverse->inverse = NULL;
880 free32(new_face);
881 return(0);
882 }
883
884 /* put a vertex in the middle of each cut edge */
885 cut_edges(polyhedron);
886
887 /* set the new face */
888 new_face->order = 0;
889 new_face->fill_tone = -1;
890 proj_copy(new_face->group_element, matrix);
891
892 /* install the new face in the polyhedron */
893 adjust_f_e_ptrs(polyhedron);
894 cut_faces(polyhedron, new_face);
895 remove_dead_edges(polyhedron);
896 remove_dead_vertices(polyhedron);
897
898 /* put the new face on the face lists */
899 new_face->next = polyhedron->face_list;
900 polyhedron->face_list = new_face;
901 new_face->nxt = polyhedron->dirty0.nxt;
902 new_face->prv = &polyhedron->dirty0;
903 polyhedron->dirty0.nxt->prv = new_face;
904 polyhedron->dirty0.nxt = new_face;
905
906 polyhedron->num_faces++;
907
908 return(1);
909 }
910
911
cut_edges(polyhedron)912 static void cut_edges(polyhedron)
913 WEpolyhedron *polyhedron;
914 {
915 int i;
916 double d0,
917 d1,
918 t,
919 s;
920 WEedge *edge,
921 *new_edge,
922 *nbr_edge;
923 WEvertex *new_vertex;
924
925 for (edge=polyhedron->edge_list; edge; edge=edge->next) {
926 d0 = edge->v0->dist;
927 d1 = edge->v1->dist;
928 if ((d0 < 0.0 && d1 > 0.0) || (d0 > 0.0 && d1 < 0.0)) {
929 new_vertex = (WEvertex *) malloc32((MyInt32) sizeof(WEvertex));
930 new_edge = (WEedge *) malloc32((MyInt32) sizeof(WEedge));
931
932 t = -d0/(d1 - d0);
933 s = 1.0 - t;
934 for (i=4; --i>=0; )
935 new_vertex->x[i] = s * edge->v0->x[i] + t * edge->v1->x[i];
936 new_vertex->dist = 0.0;
937 new_vertex->next = polyhedron->vertex_list;
938 polyhedron->vertex_list = new_vertex;
939 polyhedron->num_vertices++;
940
941 new_edge->v1 = edge->v1;
942 new_edge->v0 = new_vertex;
943 edge->v1 = new_vertex;
944
945 nbr_edge = edge->e1L;
946 new_edge->e1L = nbr_edge;
947 if (nbr_edge->e0L == edge)
948 nbr_edge->e0L = new_edge;
949 else
950 nbr_edge->e1R = new_edge;
951
952 nbr_edge = edge->e1R;
953 new_edge->e1R = nbr_edge;
954 if (nbr_edge->e0R == edge)
955 nbr_edge->e0R = new_edge;
956 else
957 nbr_edge->e1L = new_edge;
958
959 new_edge->e0L = edge;
960 new_edge->e0R = edge;
961 edge->e1L = new_edge;
962 edge->e1R = new_edge;
963
964 new_edge->fL = edge->fL;
965 new_edge->fR = edge->fR;
966
967 edge->fL->order++;
968 edge->fR->order++;
969
970 new_edge->next = polyhedron->edge_list;
971 polyhedron->edge_list = new_edge;
972 polyhedron->num_edges++;
973 }
974 }
975 return;
976 }
977
978
adjust_f_e_ptrs(polyhedron)979 static void adjust_f_e_ptrs(polyhedron)
980 WEpolyhedron *polyhedron;
981 {
982 WEedge *edge;
983
984 /* make sure each good face "sees" a good edge */
985
986 for (edge=polyhedron->edge_list; edge; edge=edge->next)
987 if (edge->v0->dist < 0.0 || edge->v1->dist < 0.0) {
988 edge->fL->some_edge = edge;
989 edge->fR->some_edge = edge;
990 }
991
992 return;
993 }
994
995
996 /* check each old face: */
997 /* (1) if a face is entirely >= 0, remove it */
998 /* (this can be checked immediately, because each */
999 /* good face sees a good edge at this point) */
1000 /* (check that removing dead faces doesn't */
1001 /* change the group) */
1002 /* (2) if a face is entirely <=0, leave it alone */
1003 /* (but make sure any 0-0 edges "see" the new face */
1004 /* and the new face sees the 0-0 edges) */
1005 /* also set order of new face */
1006 /* (3) otherwise bisect the face with a new edge, and */
1007 /* make sure the new edge "sees" the new new face */
1008 /* and the old face "sees" a valid edge */
1009
cut_faces(polyhedron,new_face)1010 static void cut_faces(polyhedron, new_face)
1011 WEpolyhedron *polyhedron;
1012 WEface *new_face;
1013 {
1014 int zero_count,
1015 count,
1016 count1 = 0,
1017 count3 = 0;
1018 double d0, d1;
1019 WEvertex *v01, *v23;
1020 WEedge *edge,
1021 *next_edge,
1022 *new_edge,
1023 *e0 = NULL, *e1 = NULL, *e2 = NULL, *e3 = NULL;
1024 WEface *f_list,
1025 *face;
1026
1027 /* to facilitate removing unneeded faces, first move all */
1028 /* the faces onto f_list, then move the good faces back */
1029 /* onto polyhedron->face_list as they are processed */
1030 f_list = polyhedron->face_list;
1031 polyhedron->face_list = NULL;
1032
1033 /* we'll count the order of the new_face as we go */
1034 new_face->order = 0;
1035
1036 while (f_list) {
1037 /* pull a face off f_list */
1038 face = f_list;
1039 f_list = f_list->next;
1040
1041 /* Is the face entirely >= 0 ? */
1042 /* Note: adjust_f_e_ptrs() has been called, */
1043 /* so all good faces see good edges. */
1044 /* Some fL and fR pointers on 0-0 edges may */
1045 /* temporarily be left dangling when face is freed. */
1046 if (face->some_edge->v0->dist >= 0
1047 && face->some_edge->v1->dist >= 0) {
1048 if (face->inverse)
1049 face->inverse->inverse = NULL;
1050 face->prv->nxt = face->nxt;
1051 face->nxt->prv = face->prv;
1052 free32(face);
1053 --polyhedron->num_faces;
1054 continue;
1055 }
1056
1057 /* cut the face if necessary */
1058
1059 /* counts number of vertices at dist 0 (two consecutive */
1060 /* vertices at dist 0 are counted as one (or none) */
1061 /* because we don't need to cut such a face) */
1062 zero_count = 0;
1063
1064 /* We'll traverse the face counterwise to find the */
1065 /* edges going in and out of the "zero vertices" (i.e. */
1066 /* the vertices at dist 0) */
1067 /* e0 = negative to zero */
1068 /* e1 = zero to positive */
1069 /* e2 = positive to zero */
1070 /* e3 = zero to negative */
1071
1072 count = 0; /* for use in finding order of new face */
1073 edge = face->some_edge;
1074 do {
1075 /* which way does the edge point? */
1076 /* edge points counterclockwise */
1077 if (edge->fL == face) {
1078 d0 = edge->v0->dist;
1079 d1 = edge->v1->dist;
1080 next_edge = edge->e1L;
1081 }
1082 else { /* edge points clockwise */
1083 d0 = edge->v1->dist;
1084 d1 = edge->v0->dist;
1085 next_edge = edge->e0R;
1086 }
1087
1088 if (d0 == 0.0) {
1089 if (d1 == 0.0) {
1090 if (edge->fL == face)
1091 edge->fR = new_face;
1092 else
1093 edge->fL = new_face;
1094 new_face->some_edge = edge;
1095 new_face->order++;
1096 break;
1097 }
1098 zero_count++;
1099 if (d1 < 0.0) {
1100 e3 = edge;
1101 count3 = count;
1102 }
1103 else { /* d1 > 0.0 */
1104 e1 = edge;
1105 count1 = count;
1106 }
1107 }
1108 else if (d1 == 0.0) {
1109 if (d0 < 0.0)
1110 e0 = edge;
1111 else /* d0 > 0.0 */
1112 e2 = edge;
1113 }
1114
1115 edge = next_edge;
1116 count++;
1117 }
1118 while (edge != face->some_edge);
1119
1120 if (zero_count == 2) { /* we need to make a cut */
1121 new_edge = &edgeheap[heapcount++];
1122 new_edge = (WEedge *) malloc32((MyInt32) sizeof(WEedge));
1123 polyhedron->num_edges++;
1124 new_edge->next = polyhedron->edge_list;
1125 polyhedron->edge_list = new_edge;
1126
1127 /* v01 = vertex between edges e0 and e1 */
1128 /* v23 = vertex between edges e2 and e3 */
1129 v01 = (e0->v0->dist == 0.0) ? e0->v0 : e0->v1;
1130 v23 = (e2->v0->dist == 0.0) ? e2->v0 : e2->v1;
1131
1132 new_edge->v0 = v01;
1133 new_edge->v1 = v23;
1134 new_edge->e0L = e0;
1135 new_edge->e0R = e1;
1136 new_edge->e1L = e3;
1137 new_edge->e1R = e2;
1138 new_edge->fL = face;
1139 new_edge->fR = new_face;
1140
1141 if (e0->v0 == v01)
1142 e0->e0R = new_edge;
1143 else
1144 e0->e1L = new_edge;
1145
1146 if (e1->v0 == v01)
1147 e1->e0L = new_edge;
1148 else
1149 e1->e1R = new_edge;
1150
1151 if (e2->v0 == v23)
1152 e2->e0R = new_edge;
1153 else
1154 e2->e1L = new_edge;
1155
1156 if (e3->v0 == v23)
1157 e3->e0L = new_edge;
1158 else
1159 e3->e1R = new_edge;
1160
1161 new_face->some_edge = new_edge;
1162 new_face->order++;
1163
1164 face->order = (count1 - count3 + face->order)%face->order + 1;
1165
1166 /* transfer face to dirty list */
1167 face->prv->nxt = face->nxt;
1168 face->nxt->prv = face->prv;
1169 face->nxt = polyhedron->dirty0.nxt;
1170 face->prv = &polyhedron->dirty0;
1171 polyhedron->dirty0.nxt->prv = face;
1172 polyhedron->dirty0.nxt = face;
1173 }
1174
1175 /* put the face back on polyhedron->face_list */
1176 face->next = polyhedron->face_list;
1177 polyhedron->face_list = face;
1178 }
1179
1180 return;
1181 }
1182
1183
remove_dead_edges(polyhedron)1184 static void remove_dead_edges(polyhedron)
1185 WEpolyhedron *polyhedron;
1186 {
1187 WEedge *e_list,
1188 *edge;
1189
1190 /* Move all edges onto e_list, then move the good edges */
1191 /* back onto polyhedron->edge_list. */
1192
1193 e_list = polyhedron->edge_list;
1194 polyhedron->edge_list = NULL;
1195
1196 while (e_list) {
1197 /* pull an edge off e_list */
1198 edge = e_list;
1199 e_list = e_list->next;
1200
1201 /* if it has one or more bad vertices, throw it away */
1202 if (edge->v0->dist > 0.0 || edge->v1->dist > 0.0) {
1203
1204 /* first tidy up at vertices of dist 0.0 */
1205 if (edge->v0->dist == 0.0) {
1206 if (edge->e0L->e0R == edge)
1207 edge->e0L->e0R = edge->e0R;
1208 else
1209 edge->e0L->e1L = edge->e0R;
1210
1211 if (edge->e0R->e0L == edge)
1212 edge->e0R->e0L = edge->e0L;
1213 else
1214 edge->e0R->e1R = edge->e0L;
1215 }
1216 if (edge->v1->dist == 0.0) {
1217 if (edge->e1L->e1R == edge)
1218 edge->e1L->e1R = edge->e1R;
1219 else
1220 edge->e1L->e0L = edge->e1R;
1221
1222 if (edge->e1R->e1L == edge)
1223 edge->e1R->e1L = edge->e1L;
1224 else
1225 edge->e1R->e0R = edge->e1L;
1226 }
1227
1228 /* throw away the edge */
1229 free32(edge);
1230 --polyhedron->num_edges;
1231 }
1232 else { /* put it back on polyhedron->edge_list */
1233 edge->next = polyhedron->edge_list;
1234 polyhedron->edge_list = edge;
1235 }
1236 }
1237 return;
1238 }
1239
1240
remove_dead_vertices(polyhedron)1241 static void remove_dead_vertices(polyhedron)
1242 WEpolyhedron *polyhedron;
1243 {
1244 WEvertex *v_list,
1245 *vertex;
1246
1247 /* Move all vertices onto v_list, then move the */
1248 /* good vertices back onto polyhedron->vertex_list. */
1249
1250 v_list = polyhedron->vertex_list;
1251 polyhedron->vertex_list = NULL;
1252
1253 while (v_list) {
1254 /* pull an vertex off v_list */
1255 vertex = v_list;
1256 v_list = v_list->next;
1257
1258 /* if it's been cut off, throw it away */
1259 if (vertex->dist > 0.0) {
1260 free32(vertex);
1261 --polyhedron->num_vertices;
1262 }
1263 else { /* put it back on polyhedron->vertex_list */
1264 vertex->next = polyhedron->vertex_list;
1265 polyhedron->vertex_list = vertex;
1266 }
1267 }
1268 return;
1269 }
1270
1271
number_faces(polyhedron)1272 static void number_faces(polyhedron)
1273 WEpolyhedron *polyhedron;
1274 {
1275 int count;
1276 WEface *face;
1277
1278 /* at this point just assign consecutive indices to the */
1279 /* faces, and let the routines for making display lists,*/
1280 /* saving in other file formats, etc. set the actual */
1281 /* fill tones */
1282
1283 count = 0;
1284 for (face=polyhedron->face_list; face; face=face->next) {
1285 /* skip faces which were already assigned fill_tones */
1286 /* when their inverses were found */
1287 if (face->fill_tone >= 0)
1288 continue;
1289
1290 if (face->inverse == NULL) {
1291 if (debug) fprintf(stderr, "unmatched faces in Dirichlet.c\n");
1292 }
1293
1294 face->fill_tone = count;
1295 if (face->inverse) face->inverse->fill_tone = count;
1296 count++;
1297 }
1298
1299 return;
1300 }
1301
1302 #if 0
1303 static void read_vertices(face, v)
1304 WEface *face;
1305 double (*v)[3];
1306 {
1307 int i;
1308 WEedge *edge;
1309
1310 edge = face->some_edge;
1311 do {
1312 if (edge->fL == face) { /* edge points counterclockwise */
1313 for (i=3; --i>=0; )
1314 (*v)[i] = edge->v0->x[i];
1315 v++;
1316 edge = edge->e1L;
1317 }
1318 else { /* edge points clockwise */
1319 for (i=3; --i>=0; )
1320 (*v)[i] = edge->v1->x[i];
1321 v++;
1322 edge = edge->e0R;
1323 }
1324 }
1325 while (edge != face->some_edge);
1326
1327 return;
1328 }
1329 #endif
1330
print_poly(polyhedron)1331 static void print_poly(polyhedron)
1332 WEpolyhedron *polyhedron;
1333 {
1334 print_vef(polyhedron);
1335 print_vertices(polyhedron);
1336 }
1337
1338 #if 0
1339 static void print_statistics(polyhedron)
1340 WEpolyhedron *polyhedron;
1341 {
1342 print_vef(polyhedron);
1343 print_vertex_distances(polyhedron);
1344 print_edge_lengths(polyhedron);
1345 print_face_distances(polyhedron);
1346 return;
1347 }
1348 #endif
1349
print_vef(polyhedron)1350 static void print_vef(polyhedron)
1351 WEpolyhedron *polyhedron;
1352 {
1353 if (debug) fprintf(stderr, "%d vertices, %d edges, %d faces\n",
1354 polyhedron->num_vertices,
1355 polyhedron->num_edges,
1356 polyhedron->num_faces);
1357
1358 if (polyhedron->num_vertices - polyhedron->num_edges + polyhedron->num_faces != 2) {
1359 if (debug) fprintf(stderr, "Euler characteristic error in Dirichlet.c\n");
1360 return;
1361 }
1362
1363 return;
1364 }
1365
1366 #if 0
1367 static void saveOOGL(polyhedron)
1368 WEpolyhedron *polyhedron;
1369 {
1370 FILE *fp = fopen("/tmp/test.off","w");
1371 GeomFSave(WEPolyhedronToPolyList(polyhedron), fp, "/tmp/test.off");
1372 fclose(fp);
1373 }
1374 #endif
1375
print_vertices(polyhedron)1376 static void print_vertices(polyhedron)
1377 WEpolyhedron *polyhedron;
1378 {
1379 WEvertex *vertex;
1380 fprintf(stderr, "Vertices:\n");
1381 for (vertex=polyhedron->vertex_list; vertex; vertex=vertex->next)
1382 fprintf(stderr, "%f\t%f\t%f\t%f\n",vertex->x[0], vertex->x[1], vertex->x[2], vertex->x[3]);
1383 }
1384
1385 #if 0
1386 static void print_vertex_distances(polyhedron)
1387 WEpolyhedron *polyhedron;
1388 {
1389 WEvertex *vertex;
1390 int ideal_vertex_present,
1391 finite_vertex_present,
1392 i;
1393 double norm,
1394 max,
1395 min,
1396 ideal_point_norm;
1397
1398 /* Note: the following statement revealed that "infinity" */
1399 /* is at a distance of about 17 on a Mac. */
1400 /* fprintf(stderr, "infinity = %lf\n", acosh(1.0/ideal_point_norm)); */
1401
1402 ideal_point_norm = sqrt(1e5 * HARDWARE_PRECISION);
1403 ideal_vertex_present = 0;
1404 finite_vertex_present = 0;
1405 min = cosh(17.0);
1406 max = 0.0;
1407
1408 if (metric == DG_HYPERBOLIC) {
1409 for (vertex=polyhedron->vertex_list; vertex; vertex=vertex->next) {
1410 norm = sqrt(fabs(DHPt3Dot3(vertex->x, vertex->x, metric)));
1411 if (norm < ideal_point_norm) {
1412 vertex->ideal = 1;
1413 ideal_vertex_present = 1;
1414 }
1415 else {
1416 vertex->ideal = 0;
1417
1418 for (i=4; --i>=0; )
1419 vertex->x[i] /= norm;
1420
1421 if (vertex->x[3] < min)
1422 min = vertex->x[3];
1423
1424 if (vertex->x[3] > max)
1425 max = vertex->x[3];
1426
1427 finite_vertex_present = 1;
1428 }
1429 }
1430
1431 if (finite_vertex_present) {
1432 fprintf(stderr, "closest vertex %.6f\n", acosh(min));
1433 if (ideal_vertex_present)
1434 fprintf(stderr, "furthest finite vertex %.6f\n", acosh(max));
1435 else
1436 fprintf(stderr, "furthest vertex %.6f\n", acosh(max));
1437 }
1438 else
1439 fprintf(stderr, "all vertices are at infinity\n");
1440
1441 }
1442 return;
1443 }
1444 #endif
1445
1446 /* Note: print_edge_lengths() assumes print_vertex_distances() */
1447 /* has already been called to normalize the vertices. */
1448
1449 #if 0
1450 static void print_edge_lengths(polyhedron)
1451 WEpolyhedron *polyhedron;
1452 {
1453 WEedge *edge;
1454 int finite_edge_present,
1455 infinite_edge_present;
1456 double dot,
1457 min,
1458 max;
1459
1460 min = cosh(17.0);
1461 max = 1.0;
1462 finite_edge_present = 0;
1463 infinite_edge_present = 0;
1464
1465 for (edge=polyhedron->edge_list; edge; edge=edge->next)
1466 if (edge->v0->ideal || edge->v1->ideal)
1467 infinite_edge_present = 1;
1468 else {
1469 finite_edge_present = 1;
1470 dot = -DHPt3Dot3( edge->v0->x,edge->v1->x, metric);
1471 if (dot < min)
1472 min = dot;
1473 if (dot > max)
1474 max = dot;
1475 }
1476
1477 if (finite_edge_present)
1478 if (infinite_edge_present) {
1479 fprintf(stderr, "shortest finite edge %.6f\n", acosh(min));
1480 fprintf(stderr, "longest finite edge %.6f\n", acosh(max));
1481 }
1482 else {
1483 fprintf(stderr, "shortest edge %.6f\n", acosh(min));
1484 fprintf(stderr, "longest edge %.6f\n", acosh(max));
1485 }
1486 else
1487 fprintf(stderr, "all edges are infinite\n");
1488
1489 return;
1490 }
1491 #endif
1492
1493 #if 0
1494 static void print_face_distances(polyhedron)
1495 WEpolyhedron *polyhedron;
1496 {
1497 WEface *face;
1498 double min,
1499 max;
1500
1501 min = 1e308; /* ieee max double */
1502 max = 0;
1503
1504 for (face=polyhedron->face_list; face; face=face->next) {
1505 if (face->group_element[3][3] < min)
1506 min = face->group_element[3][3];
1507 if (face->group_element[3][3] > max)
1508 max = face->group_element[3][3];
1509 }
1510
1511 fprintf(stderr, "closest face plane %.6f\n", 0.5 * acosh(min));
1512 fprintf(stderr, "furthest face plane %.6f\n", 0.5 * acosh(max));
1513
1514 return;
1515 }
1516 #endif
1517
1518 #if 0
1519 static void free_polyhedron(polyhedron)
1520 WEpolyhedron *polyhedron;
1521 {
1522 WEvertex *dead_vertex;
1523 WEedge *dead_edge;
1524 WEface *dead_face;
1525
1526 while (polyhedron->vertex_list) {
1527 dead_vertex = polyhedron->vertex_list;
1528 polyhedron->vertex_list = dead_vertex->next;
1529 free32(dead_vertex);
1530 }
1531
1532 while (polyhedron->edge_list) {
1533 dead_edge = polyhedron->edge_list;
1534 polyhedron->edge_list = dead_edge->next;
1535 free32(dead_edge);
1536 }
1537
1538 while (polyhedron->face_list) {
1539 dead_face = polyhedron->face_list;
1540 polyhedron->face_list = dead_face->next;
1541 free32(dead_face);
1542 }
1543
1544 return;
1545 }
1546 #endif
1547
roundoff_message(epsilon_name)1548 static void roundoff_message(epsilon_name)
1549 char *epsilon_name;
1550 {
1551 #if PRECISE_GENERATORS
1552 fprintf(stderr,"\nWARNING: roundoff error is getting a bit large. (%s)\n", epsilon_name);
1553 #else
1554 fprintf(stderr,"\nWARNING: roundoff error is getting perilously large. (%s)\n", epsilon_name);
1555 #endif
1556 fprintf(stderr,"To verify the correctness of the final Dirichlet domain,\n");
1557 fprintf(stderr,"move off to another point in the Dehn filling plane, then\n");
1558 fprintf(stderr,"return to the present point, recompute the Dirichlet domain,\n");
1559 fprintf(stderr,"and see whether you get the same number of vertices, edges\n");
1560 fprintf(stderr,"and faces. (Moving to a different point and then returning\n");
1561 fprintf(stderr,"will randomize the roundoff error.)\n\n");
1562
1563 return;
1564 }
1565
proj_same_matrix(m0,m1)1566 int proj_same_matrix(m0, m1)
1567 proj_matrix m0,
1568 m1;
1569 {
1570 int i, j;
1571 double diff;
1572
1573 for (i=4; --i>=0; )
1574 for (j=4; --j>=0; ) {
1575 diff = fabs(m0[i][j] - m1[i][j]);
1576 if (diff > MATRIX_EPSILON)
1577 return(0);
1578 if (diff > 1e-2 * MATRIX_EPSILON && ! matrix_epsilon_message_given) {
1579 if (debug) roundoff_message("MATRIX_EPSILON");
1580 matrix_epsilon_message_given = 1;
1581 }
1582 }
1583 return(1);
1584 }
1585