conformalmesher.C

Go to the documentation of this file.

#ifdef HAVE_CONFIG_H
#include <tumble-conf.h>
#endif /* HAVE_CONFIG_H */

#include "conformalmesher.h"

#include "beziermesh.h"
#include "boundarymesh.h"
#include "cell.h"
#include "globals.h"
#include "iterator_wrapper.h"
#include "spline.h"

#include <algorithm>
#include <stack>


/* Turn on verbose messages from the conformal mesher */
//#define DEBUG_CONF_MESHER

#ifdef DEBUG_CONF_MESHER
  #include "meshio.h"
#endif

using namespace std;
using namespace hashers;

class ConformalMesher::Subspline {
  public:
    Subspline(BoundaryEdge *bdrye, BezierVertex *v0, BezierVertex *v1) :
      bdrye_(bdrye) {
        v_[0] = v0;
        v_[1] = v1;
      }

    /* for O(1) destruction without leaking memory. */
    void set_delete_iterator(std::list<Subspline*>::iterator it) {
      delete_iterator = it;
    }
    list<Subspline*>::iterator get_delete_iterator() {
      return delete_iterator;
    }

    BezierVertex *get_vertex(int i) {
      assert(i >= 0); assert(i < 2);
      return v_[i];
    }
    const BezierVertex *get_vertex(int i) const {
      return const_cast<Subspline*>(this)->get_vertex(i);
    }

    const BezierVertex *get_other_vertex(BezierVertex *v) const {
      if(get_vertex(0) == v) {
        return get_vertex(1);
      } else {
        assert(get_vertex(1) == v);
        return get_vertex(0);
      }
    }

    const QBSpline *get_spline() const { return bdrye_->get_spline(); }
    QBSpline *get_spline() { return bdrye_->get_spline(); }

    const Point2D& get_deboor() const {
      int idx = get_segment_index();
      return *get_spline()->get_deboor(idx + 1); /* segment 0 has deboor 1 */
    }

    BezierVertex *split(BezierMesh *mesh, vector<Subspline*>& subs) {
      /* split the spline; return two new subsplines.
         the caller is responsible for deleting this subspline */
      int idx = get_segment_index();
      ControlPoint cp[3];
      ControlPoint newpt;

      get_spline()->split(idx, 0.5, cp[0], cp[1], cp[2]);
      newpt = cp[1]; /* we ignore cp[0] and cp[2] */

      BezierVertex *splitvertex = mesh->clean_insert_point(*newpt);
      splitvertex->set_bdry(bdrye_, idx + 1);
      mesh->replace_control_point(splitvertex, newpt);

      /* order of the vertices is important */
      subs.push_back(new Subspline(bdrye_, get_vertex(0), splitvertex));
      subs.push_back(new Subspline(bdrye_, splitvertex, get_vertex(1)));

      return splitvertex;
    }

    /* check if the vertex encroaches on us */
    bool is_encroached(BezierVertex *v) const {
      int idx = get_segment_index();
      /* '2' here is the degree of the curve: segment 0 starts
       * at b[0] and goes until b[2]. */
      Tumble::iterators::dereferencer<typeof(get_spline()->b.begin()), Point2D,
    const Point2D&>
       first(get_spline()->b.begin() + 2 * idx),
       last (get_spline()->b.begin() + 2 * (idx + 1));
      return BezierEdge::encroaches(v->get_cp(), first, ++last);
    }

    int get_segment_index() const {
      if(get_vertex(0)->containment_dimension() == 0) {
        /* our left vertex is a D0 vertex? then we're the first segment */
        return 0;
      }
      else {
        int idx = get_spline()->get_bezier_vertex_idx(get_vertex(0));
        assert(idx > 0); assert(idx < get_spline()->get_num_segments());
        return idx;
      }
    }


  private:
    BoundaryEdge *bdrye_;
    BezierVertex *v_[2];
    list<Subspline*>::iterator delete_iterator;
};

#ifdef DEBUG_CONF_MESHER
static ostream& operator<<(ostream& os, ConformalMesher::Subspline *s) {
  return os << (void*)s << " [spline: " << (s->get_spline()) << " => "
    << s->get_vertex(0) << " to " << s->get_vertex(1)
    << " via " << s->get_deboor() << ']';
}
#endif

/***************************************************************************
  ...
 ***************************************************************************/
ConformalMesher::ConformalMesher(BezierMesh* _bezier_mesh,
                                 BoundaryMesh* _bdry_mesh)
{
  bezier = _bezier_mesh;
  bdry = _bdry_mesh;
}

/***************************************************************************
  Delete all the subsplines.
 ***************************************************************************/
ConformalMesher::~ConformalMesher()
{
  for(list<Subspline*>::iterator it = subsplines.begin();
      it != subsplines.end(); ++it) {
    delete (*it);
  }
}


/***************************************************************************
  The main method to build a conformal mesh.  Simply passes off the work
  to others.
 ***************************************************************************/
void ConformalMesher::mesh()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: starting.\n";
#endif
  /* everything is linear at the start; allow the mesh to assume so */
  bezier->set_linear_flips(true);

  create_bbox();
  insert_d0_vertices();
  create_subsplines();
  conform();
  reticulate_splines();

  /* The mesh now is fully built; we're down to painting. */
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: done building; mesh size "
    << bezier->get_num_faces() << endl;
  cout<<"BezierMesh:"<<endl;
  bezier->print_statistics();
#endif
  bezier->set_linear_flips(false);

  set_face_pointers();
  excavate_holes();
  cout << "ConformalMesher: done.\n";
  cout<<"************ Statistics (make sure these add up) ****************"<<endl;
  cout<<"BoundaryMesh:"<<endl;
  bdry->print_statistics();
  cout<<"BezierMesh:"<<endl;
  bezier->print_statistics();
}

/***************************************************************************
  Create a bounding box and triangulate it.
 ***************************************************************************/
void ConformalMesher::create_bbox()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: building bbox" << flush;
#endif
  /* [x][y] -- so [0][1] is the top-left, [1][0] bottom-right */
  Point2D p[2][2];
  BezierVertex *v[2][2];
  BezierEdge *e[5]; /* 4 outside, 1 cross edge */
  double diameter;

  bdry->bbox(p[1][1], p[0][0]);
  p[1][0] = Point2D(p[1][1].x(), p[0][0].y());
  p[0][1] = Point2D(p[0][0].x(), p[1][1].y());

  /* Grow the v such that there is no way that the exterior is
     introducing a small feature. */
  diameter = (p[0][0]-p[1][1]).mag();
  p[0][0] += Point2D(-diameter, -diameter);
  p[0][1] += Point2D(-diameter,  diameter);
  p[1][0] += Point2D( diameter, -diameter);
  p[1][1] += Point2D( diameter,  diameter);
  cout << " bottom-left: " << p[0][0] << " top-right " << p[1][1] << endl;

  /* Create the vertices. */
  for(int i = 0; i < 2; i++) {
    for(int j = 0; j < 2; j++) {
      v[i][j] = bezier->add_bezier_vertex(p[i][j]);
    }
  }

  /* Create the exterior and the one cross edge. */
  e[0] = bezier->add_straight_bezier_edge(v[0][0], v[0][1]);
  e[1] = bezier->add_straight_bezier_edge(v[0][0], v[1][0]);
  e[2] = bezier->add_straight_bezier_edge(v[1][1], v[0][1]);
  e[3] = bezier->add_straight_bezier_edge(v[1][1], v[1][0]);
  e[4] = bezier->add_straight_bezier_edge(v[0][0], v[1][1]);

  /* Create the two faces, in ccw order. */
  bezier->add_bezier_triangle(NULL, e[4], e[2], e[0]); //CCW: tp517 4/16/07
  bezier->add_bezier_triangle(NULL, e[4], e[1], e[3], true); // CCW: tp517 4/16/07

  /* The mesh is now Delaunay.  Remember the need to excavate the
   * exterior face when we're done. */
  remove_points.push_back(p[0][0]);

#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}

/***************************************************************************
  Insert the D0 vertices (input vertices).
  Just vanilla insertions; should be random but we the meshes are assumed
  to be small.
 ***************************************************************************/
void ConformalMesher::insert_d0_vertices()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: inserting D0 vertices\n";
#endif
  for(BoundaryMesh::Boundary_Vertex_Hash_T::iterator it = bdry->get_vertexs_begin();
      it != bdry->get_vertexs_end(); ++it) {
    BoundaryVertex *bdryv = (*it);
    BezierVertex *bezv = bezier->clean_insert_point(bdryv->get_cp());
    bezv->set_bdry(bdryv);
    bezier->replace_control_point(bezv, bdryv->get_control_point());

#ifdef DEBUG_CONF_MESHER
    cout << "ConformalMesher:     - inserted " << (void*)bezv
      << " at " << bezv->get_cp() << endl;
#endif
  }

#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}

/***************************************************************************
  Create the splines:
  - insert D1 vertices (the knots)
  - create Subsplines and put them on the stack
 ***************************************************************************/
void ConformalMesher::create_subsplines()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: inserting D1 vertices, storing splines\n";
#endif
  for(BoundaryMesh::Boundary_Edge_Hash_T::iterator it = bdry->get_edges_begin();
      it != bdry->get_edges_end(); ++it) {
    BoundaryEdge *bdrye = (*it);
    QBSpline *spline = bdrye->get_spline();
    vector<BezierVertex*> knots;

#ifdef DEBUG_CONF_MESHER
    vector<BezierTriangle*> inverts;
#endif

    /* Incrementally add every d1 vertex of the spline; store
     * all the d0 and d1 vertices, correctly indexed, in 'knots'.
     */
    knots.push_back(bdrye->get_vertex(0)->get_bezier());
    for (int i = 1; i <= spline->get_num_segments() - 1; i++) {
      ControlPoint cp = spline->get_control_point_at_idx(i);
      BezierVertex *bezv = bezier->clean_insert_point(*cp);

#ifdef DEBUG_CONF_MESHER
      cout << "ConformalMesher:     - inserting " << cp << "... " << flush;
      bezier->find_inverted_triangles(inverts);
      if (!inverts.empty())
    {
      cout<< " INVERTED TRIANGLES FOUND IN CONF_MESHER "<<endl;
      cout<< " Just inserted : "<<cp<<endl;
      for (vector<BezierTriangle*>::iterator i = inverts.begin();
           i != inverts.end(); ++i)
        cout << "Inverted Triangle: "<<(**i)<<endl;
      MeshBinaryOutput mo(bezier,bdry);
      mo.write("created_inverted_during_create_subsplines.geo");
      assert(0);
    } else {
      cout << "Conformal mesher inversion check passed."<<endl;
    }
#endif

      bezier->replace_control_point(bezv, cp);
      bezv->set_bdry(bdrye, i);
      knots.push_back(bezv);

      spline->set_bezier_vertex(i, knots.back());

    }
    knots.push_back(bdrye->get_vertex(1)->get_bezier());
    assert((unsigned)knots.size() == (unsigned)spline->get_num_segments() + 1);

    /* Create every subspline but don't put them in the mesh yet. */
    for (int i = 0; i < spline->get_num_segments(); i++) {
      Subspline *s = new Subspline(bdrye, knots[i], knots[i+1]);
      link_edge(s);
    }
  }
#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}


/***************************************************************************
  Enforce conformity:
  while there are subsplines that either don't appear or are encroached,
      split the subspline with a new vertex
 ***************************************************************************/
void ConformalMesher::conform()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: enforcing conformity\n";
#endif
  while(!spline_stack.empty()) {
    Subspline *s = spline_stack.back();
    spline_stack.pop_back();
    in_stack.erase(s);
#ifdef DEBUG_CONF_MESHER
    cout << "ConformalMesher:   - considering " << s << endl;
    cout << "ConformalMesher:    ";
#endif
    if ( (is_encroached.count(s) == 0) && segment_appears(s)) {
      #ifdef DEBUG_CONF_MESHER
      cout << "+ happy\n";
      #endif
      continue;
    }
#ifdef DEBUG_CONF_MESHER
    cout << " => splitting!\n";
#endif

    vector<Subspline*> newsegs;
    BezierVertex *splitvertex = split_destroy_segment(s, newsegs);
    assert(newsegs.size() == 2);
    s = (Subspline*)0x7; /* for debugging, make it crash */
#ifdef DEBUG_CONF_MESHER
    cout << "ConformalMesher:      produced:\n"
         << "                      " << newsegs[0] << endl
         << "                      " << newsegs[1] << endl;
#endif

    /* check neighbours of the new vertex for encroachment */
    deque<BezierEdge*> link_edges;
    bezier->enqueue_edges(splitvertex, link_edges);
    for(deque<BezierEdge*>::iterator link_it = link_edges.begin();
        link_it != link_edges.end(); ++link_it) {
      BezierEdge *link_e = (*link_it);
      BezierVertex *link_v = bezier->get_opposite_vertex(link_e, splitvertex);

      vector<Subspline*>::iterator seg_it = incident_edges[link_v].begin();
      vector<Subspline*>::iterator seg_end = incident_edges[link_v].end();
      for ( ; seg_it != seg_end ; ++seg_it) {
        Subspline *link_s = *seg_it;
        if (link_s == newsegs[0] || link_s == newsegs[1]) {
          /* splitvertex doesn't encroach on the segments of which it's
             a member */
          continue;
        }
        if (is_encroached.count(link_s) != 0) {
          /* Already both in the stack, and marked.
             Don't quit if !is_encroached but in_stack: when we pop the
             segment off, we must know about the encroachment */
          assert(in_stack.count(link_s) != 0);
          continue;
        }

        /* Check for encroachment; mark and ensure it's on the stack if it
         * is encroached. */
        if (link_s->is_encroached(splitvertex)) {
          is_encroached.insert(link_s);
          if (in_stack.count(link_s) == 0) {
            spline_stack.push_back(link_s);
            in_stack.insert(link_s);
          }
        }
      }
    }
  }

#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}

/***************************************************************************
  Seek out and 'fix' the BezierEdges that think they're straight but are
  actually conforming to a boundary.  Assumes that conform() has been called.
 ***************************************************************************/
void ConformalMesher::reticulate_splines()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: adding splines to BezierMesh\n";
#endif

  for (BoundaryMesh::Boundary_Edge_Hash_T::iterator it = bdry->get_edges_begin();
      it != bdry->get_edges_end(); ++it) {
    BoundaryEdge *bdrye = (*it);
    QBSpline *spline = bdrye->get_spline();

    int n = spline->get_num_segments();
    for (int i = 0 ; i < n; i++) {
      BezierVertex *v1 = (i == 0) ? bdrye->get_vertex(0)->get_bezier()
        : spline->get_bezier_vertex(i);
      assert(v1);
      BezierVertex *v2 = (i + 1 == n) ? bdrye->get_vertex(1)->get_bezier()
        : spline->get_bezier_vertex(i + 1);
      assert(v2);
      BezierEdge *edge = bezier->find_common_edge(v1, v2);
      assert(edge);

      edge->set_bdry(bdrye, spline->k[i], spline->k[i+1]);
      bezier->replace_control_point(edge, spline->get_deboor(i+1));
    }
  }
#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}

/***************************************************************************
  For every BezierTriangle corresponding to the BoundaryFace, set the
  BoundaryFace pointer.  This is where we use the face_points.
 ***************************************************************************/
void ConformalMesher::set_face_pointers()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: setting face pointers\n";
#endif
  for(vector<pair<BoundaryFace*,Point2D> >::iterator it = face_points.begin();
      it != face_points.end(); ++it) {
    BoundaryFace *face = (*it).first;
    Point2D p = (*it).second;
    hash_set<BezierTriangle*> triangles;
    int triangles_set = 0;

    collect_faces(p, triangles);
    for(hash_set<BezierTriangle*>::iterator it = triangles.begin();
        it != triangles.end(); ++it, triangles_set++) {
      BezierTriangle *t = *it;
      t->set_boundary_face(face);
    }
#ifdef DEBUG_CONF_MESHER
    cout << "ConformalMesher:   - face " << (void*)face << " set "
      << triangles_set << " triangles\n";
#endif
  }
}

/***************************************************************************
  Remove the elements in the holes.
 ***************************************************************************/
void ConformalMesher::excavate_holes()
{
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher: excavating holes\n";
#endif
  /* A list of things to kill; we must kill them in order.
     We may as well collect all holes together at once; there's
     no cost to doing otherwise.  */
  hash_set<BezierTriangle*> triangles;
  hash_set<BezierEdge*> edges;
  hash_set<BezierVertex*> vertices;

#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher:   - considering " << remove_points.size()
    << " holes\n";
#endif

  for(vector<Point2D>::iterator it = remove_points.begin();
      it != remove_points.end(); ++it) {
    collect_faces(*it, triangles);
  }

  for(hash_set<BezierTriangle*>::iterator it = triangles.begin();
      it != triangles.end(); ++it) {
    BezierTriangle *t = (*it);

    /* Gather the three edges and vertices. */
    BezierTuple final_tup = bezier->get_tuple(t);
    BezierTuple tup = final_tup;
    do {
      /* Things that aren't boundaries get removed. */
      if (!tup.e->is_boundary()) {
        edges.insert(tup.e);
      }
      if (!tup.v->is_boundary()) {
        vertices.insert(tup.v);
      }
      tup = bezier->Switch(0, bezier->Switch(1, tup));
    } while (tup != final_tup);
  }


  /* remove all triangles */
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher:   - removing " << triangles.size()
    << " triangles\n";
#endif
  for(hash_set<BezierTriangle*>::iterator it = triangles.begin();
      it != triangles.end(); ++it) {
    bezier->delete_face(*it);
  }


  /* all edges */
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher:   - removing " << edges.size() << " edges\n";
#endif
  for(hash_set<BezierEdge*>::iterator it = edges.begin();
      it != edges.end(); ++it) {
    bezier->delete_edge(*it);
  }


  /* and all vertices */
#ifdef DEBUG_CONF_MESHER
  cout << "ConformalMesher:   - removing " << vertices.size() << " vertices\n";
#endif
  for(hash_set<BezierVertex*>::iterator
      it = vertices.begin(); it != vertices.end(); ++it) {
    bezier->delete_vertex(*it);
  }

#ifdef DEBUG_CONF_MESHER
  bezier->check_consistency();
#endif
}

/***************************************************************************
 * Return whether the given segment appears in the linear mesh.
 * To appear, two things must be true:
 * (1) it appears topologically
 * (2) the curve does not intersect the sides (except at endpoints)
 ***************************************************************************/
bool ConformalMesher::segment_appears(Subspline *seg)
{
  BezierEdge *edge = bezier->find_common_edge(seg->get_vertex(0),
                                              seg->get_vertex(1));
  if (!edge) {
    /* no such edge in the mesh? clearly it doesn't appear */
    cout << "(no such edge)" << flush;
    return false;
  }

  /* There are four tangent angles to check now: one for each
   * face/edge/vertex cell that correspond to our edge.
   * Due to input restrictions, two segments intersect with at least
   * 90 degrees between them, and the control angle is at least 90
   * degrees, so curved subsegments cannot intersect.  This means we
   * can conservatively assume the four other sides of the two
   * triangles are straight; only our edge is curved.
   */
  BezierTuple tuple;
  BezierEdge *e[4];
  BezierVertex *v[4];

  /* Get the two edges out of one vertex. */
  tuple = bezier->get_tuple(edge);
  v[0] = v[1] = tuple.v;
  assert(v[0] == seg->get_vertex(0) || v[0] == seg->get_vertex(1));
  e[0] = bezier->Switch(1, tuple).e;
  e[1] = bezier->Switch(1, bezier->Switch(2, tuple)).e;

  /* And out of the other. */
  tuple = bezier->Switch(0, tuple);
  v[2] = v[3] = tuple.v;
  assert(v[2] == seg->get_vertex(0) || v[2] == seg->get_vertex(1));
  assert(v[2] != v[0]);
  e[2] = bezier->Switch(1, tuple).e;
  e[3] = bezier->Switch(1, bezier->Switch(2, tuple)).e;

  for (int i = 0; i < 4; i++) {
    assert(e[i] != edge);
    if (!angle_is_positive(seg, v[i], e[i])) {
      #ifdef DEBUG_CONF_MESHER
      cout << "(angle " << i << " is negative)" << flush;
      #endif
      return false;
    }
  }

  /* all our tests came out good; the edge appears. */
  return true;
}

/***************************************************************************
 * Given a quadratic bezier segment and a straight segment that share a
 * vertex, return true if the curve is good (that is, does not
 * intersect the triangle except at endpoints).
 * This happens iff the angle at the vertex is positive.
 **************************************************************************/
bool ConformalMesher::angle_is_positive(Subspline *s,
                                        BezierVertex *shared,
                                        BezierEdge *e)
{
  Point2D pshared = shared->get_cp();
  Point2D pedge = bezier->get_opposite_vertex(e, shared)->get_cp();
  Point2D pdeboor = s->get_deboor();
  Point2D pspline = s->get_other_vertex(shared)->get_cp();

  /* if the spline is straight, the angle is good. */
  double sp_sh_d  = pdeboor.line_side_test(pspline, pshared);
  if (sp_sh_d == 0.0) {
    return true;
  }

  /* if the deboor point is on the opposite side of the spline edge
     from pedge, then the angle is good */
  double sp_sh_e = pedge.line_side_test(pspline, pshared);
  assert(sp_sh_e != 0.0); /* that would be a flat triangle */
  if (sp_sh_d * sp_sh_e < 0.0) {
    return true;
  }

  /* if the deboor point is on the same side of the straight edge
     as pspline, then the angle is good */
  double e_sh_d = pdeboor.line_side_test(pedge, pshared);
  double e_sh_sp = pspline.line_side_test(pedge, pshared);
  assert(e_sh_sp != 0.0); /* flat triangle */
  if (e_sh_d * e_sh_sp > 0.0) {
    return true;
  }

  /* anything else means the angle is bad */
  return false;
}

/***************************************************************************
  Link in a new subspline.
 ***************************************************************************/
void ConformalMesher::link_edge(Subspline *s)
{
  spline_stack.push_back(s);
  in_stack.insert(s);
  incident_edges[s->get_vertex(0)].push_back(s);
  incident_edges[s->get_vertex(1)].push_back(s);
  subsplines.push_front(s);
  s->set_delete_iterator(subsplines.begin());
}

/***************************************************************************
  Unlink and delete a subspline.  It must not be in the stack.
 ***************************************************************************/
void ConformalMesher::unlink_edge(Subspline *s)
{
  assert(in_stack.count(s) == 0);
  is_encroached.erase(s);
  for (int vi = 0; vi < 2; ++vi) {
    vector<Subspline*>& vect = incident_edges[s->get_vertex(vi)];
    bool did_find = false;
    for (unsigned i = 0; i < vect.size(); ++i) {
      if(vect[i] == s) {
        vect[i] = vect.back();
        vect.pop_back();
        did_find = true;
        break;
      }
    }
    assert(did_find);
  }
  subsplines.erase(s->get_delete_iterator());
  delete s;
}

/***************************************************************************
  Split s, destroying it but creating two new segments and
  a new vertex.
 ***************************************************************************/
BezierVertex *ConformalMesher::split_destroy_segment(Subspline *s,
    vector<Subspline*>& newsegs)
{
  /* do the split */
  BezierVertex *splitvertex = s->split(bezier, newsegs);
  assert(newsegs.size() == 2);

  /* remove s */
  unlink_edge(s);

  /* link in the new segments */
  link_edge(newsegs[0]);
  link_edge(newsegs[1]);

  return splitvertex;
}

/***************************************************************************
  Collect all triangles that are entirely inside the boundary face in which
  the point lies.
 ***************************************************************************/
void ConformalMesher::collect_faces(Point2D start,
                                    hash_set<BezierTriangle*>& faces)
{
  int cd;
  BezierTuple searchtup;
  double u, v;

  /* Find starting point for the DFS. */
  cd = bezier->blind_locate_point(start, searchtup, u, v);

  /* do a DFS over all neighboring triangles stopping at boundaries */
  deque<BezierTriangle*> todo;

  switch(cd) {
    case 0:
      /* All faces around this vertex should be included. */
      bezier->enqueue_faces(searchtup.v, todo);
      break;

    case 1:
      /* Both faces around this edge should be included. */
      todo.push_back(searchtup.f);
      todo.push_back(bezier->Switch(2, searchtup).f);
      break;

    case 2:
      todo.push_back(searchtup.f);
      break;

    default:
      assert(0);
      break;
  }

  while(!todo.empty()){
    BezierTriangle *tri = todo.back();
    todo.pop_back();

    faces.insert(tri);

    /* loop around the three edges; flip across it to get the other
       triangle if it's not a boundary */
    BezierTuple start = bezier->get_tuple(tri);
    BezierTuple tup = start;
    do {
      if(!tup.e->is_boundary()) {
        /* cross the edge and push that face on the stack */
        BezierTriangle *t = bezier->Switch(2, tup).f;
        if(faces.count(t) == 0) {
          todo.push_back(t);
        }
      }
      /* next edge along the triangle */
      tup = bezier->Switch(0, bezier->Switch(1, tup));
    } while(tup != start);
  }
}

---