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