C++ Cell_handle类代码示例

void Foam::cellShapeControlMesh::barycentricCoords
(
    const Foam::point& pt,
    barycentric& bary,
    Cell_handle& ch
) const
{
    // Use the previous cell handle as a hint on where to start searching
    // Giving a hint causes strange errors...
    ch = locate(toPoint(pt));

    if (dimension() > 2 && !is_infinite(ch))
    {
        oldCellHandle_ = ch;

        tetPointRef tet
        (
            topoint(ch->vertex(0)->point()),
            topoint(ch->vertex(1)->point()),
            topoint(ch->vertex(2)->point()),
            topoint(ch->vertex(3)->point())
        );

        bary = tet.pointToBarycentric(pt);
    }
}

  double operator() (const AdvancingFront& adv, Cell_handle& c,
                     const int& index) const
  {
    // bound == 0 is better than bound < infinity
    // as it avoids the distance computations
    if(bound == 0){
      return adv.smallest_radius_delaunay_sphere (c, index);
    }

    // If perimeter > bound, return infinity so that facet is not used
    double d  = 0;
    d = sqrt(squared_distance(c->vertex((index+1)%4)->point(),
                              c->vertex((index+2)%4)->point()));
    if(d>bound) return adv.infinity();
    d += sqrt(squared_distance(c->vertex((index+2)%4)->point(),
                               c->vertex((index+3)%4)->point()));
    if(d>bound) return adv.infinity();
    d += sqrt(squared_distance(c->vertex((index+1)%4)->point(),
                               c->vertex((index+3)%4)->point()));
    if(d>bound) return adv.infinity();

    // Otherwise, return usual priority value: smallest radius of
    // delaunay sphere
    return adv.smallest_radius_delaunay_sphere (c, index);
  }

// Returns whether the facet facet is attached to the cell
// that is used to represent facet
bool edge_attached_to(const DT_3& dt, const Edge& edge, const Facet& facet) {
  
  if(dt.is_infinite(facet)) {
    return false;
  }

  Vertex_handle v1 = edge.first->vertex(edge.second);
  Vertex_handle v2 = edge.first->vertex(edge.third);

  Cell_handle cell = facet.first;
  int i1 = facet.second;
  int i2 = cell->index(v1);
  int i3 = cell->index(v2);
  CGAL_assertion(i1!=i2);
  CGAL_assertion(i1!=i3);
  CGAL_assertion(i2!=i3);
  int j = 0;
  
  while(j==i1 || j==i2 || j==i3) {
    j++;
  }
  // j is the index of the third point of the facet
  Vertex_handle w = cell->vertex(j);
  
  return CGAL::side_of_bounded_sphere(v1->point(),v2->point(),w->point())==CGAL::ON_BOUNDED_SIDE;

}

void draw_facets(Delaunay &Tr,std::vector<Facet> &facets,PointColor pcolors,CGAL::Geomview_stream &gv) {
  if(! gv_on) return;
  CGAL::Color colors[] = {CGAL::BLUE,CGAL::GREEN,CGAL::YELLOW,CGAL::DEEPBLUE,
			  CGAL::PURPLE,CGAL::VIOLET,CGAL::ORANGE,CGAL::RED};
  if(pcolors.size() == 0) 
    return draw_facets(Tr,facets,gv);
  // draw with color interpolation
  for (std::vector<Facet>::iterator it = facets.begin();it != facets.end();it++) {
    Facet f = *it;
    Cell_handle c = f.first;
    int j = f.second;
    CGAL::Color vcolors[3];
    int k = 0;
    for(int i = 0;i < 4;i++) 
      if(i != j)
	vcolors[k++] = pcolors[c->vertex(i)->info()];
    int r = (vcolors[0].red() + vcolors[1].red() + vcolors[2].red()) / 3;
    int g = (vcolors[0].green() + vcolors[1].green() + vcolors[2].green()) / 3;
    int b = (vcolors[0].blue() + vcolors[1].blue() + vcolors[2].blue()) / 3;
    gv << CGAL::Color(r,g,b);
    //    std::cout << "RGB " << r << " " << g << " " << b <<std::endl;
    Triangle t = Tr.triangle(f);
    gv << t;
  }
}

AlphaSimplex3D::
AlphaSimplex3D(const Delaunay3D::Edge& e)
{
    Cell_handle c = e.first;
    Parent::add(c->vertex(e.second));
    Parent::add(c->vertex(e.third));
}

double
cell_volume(const Cell_handle& c)
{
   Tetrahedron t = Tetrahedron(c->vertex(0)->point(),
                               c->vertex(1)->point(),
                               c->vertex(2)->point(),
                               c->vertex(3)->point());
   return ( CGAL::to_double(CGAL::abs(t.volume()) ) );
}

Point
circumcenter(const Facet& f)
{
    Cell_handle c = f.first;
    int id = f.second;
    Point p[3];
    for(int i = 0; i < 3; i ++)
        p[i] = c->vertex((id + (i+1))%4)->point();
    return cc_tr_3(p[0], p[1], p[2]);
}

void Foam::cellShapeControlMesh::writeTriangulation()
{
    OFstream str
    (
        "refinementTriangulation_"
      + name(Pstream::myProcNo())
      + ".obj"
    );

    label count = 0;

    Info<< "Write refinementTriangulation" << endl;

    for
    (
        CellSizeDelaunay::Finite_edges_iterator e = finite_edges_begin();
        e != finite_edges_end();
        ++e
    )
    {
        Cell_handle c = e->first;
        Vertex_handle vA = c->vertex(e->second);
        Vertex_handle vB = c->vertex(e->third);

        // Don't write far edges
        if (vA->farPoint() || vB->farPoint())
        {
            continue;
        }

        // Don't write unowned edges
        if (vA->referred() && vB->referred())
        {
            continue;
        }

        pointFromPoint p1 = topoint(vA->point());
        pointFromPoint p2 = topoint(vB->point());

        meshTools::writeOBJ(str, p1, p2, count);
    }

    if (is_valid())
    {
        Info<< "    Triangulation is valid" << endl;
    }
    else
    {
        FatalErrorIn
        (
            "Foam::triangulatedMesh::writeRefinementTriangulation()"
        )   << "Triangulation is not valid"
            << abort(FatalError);
    }
}

static void printCell(const Cell_handle &cell, const Vertex_handle &infinity)
{
	DEBUG_START;
	unsigned infinityIndex = cell->index(infinity);
	for (unsigned i = 0; i < NUM_CELL_VERTICES; ++i)
		if (i != infinityIndex)
			std::cout << cell->vertex(i)->info()
				<< " ";
	std::cout << std::endl;
	DEBUG_END;
}

/**
 * Gets the material for the given point
 * @param p The point to get the material for
 * @return The material at the given point
 */
Contents StoneWeatherer::getContents( const Point & p ) const {

	// Get a handle to the cell that contains the given point
	Cell_handle ch = newDT->locate( p );

	if ( newDT->is_infinite( ch ) ) {

		return AIR;
	}
	else {

		return ch->info();
	}
}

// Calculates void volume and spheres area (only for own four spheres) of regular triangulation cell.
// Doesn't consider extraneous atoms from neighbor cells! True volume can be obtained only for cluster of cells comprising one connected void.
// S. Sastry, D.S. Corti, P.G. Debenedetti, F.H. Stillinger, "Statistical geometry of particle packings. I. Algorithm for exact determination of connectivity,
// volume, and surface areas of void space in monodisperse and polydisperse sphere packings", Phys. Rev. E, V.56, N5, p. 5524, 1997. doi:10.1103/PhysRevE.56.5524
double cell_void_volume(Cell_handle c, double &out_surf, Array_double_4 &out_atom_surf)
{
    const Point *pts[4] = { &((c->vertex(0)->point()).point()),
                            &((c->vertex(1)->point()).point()),
                            &((c->vertex(2)->point()).point()),
                            &((c->vertex(3)->point()).point()) };
    const Point &V = c->weighted_circumcenter().point();   // Voronoi vertex V

    double volume = 0.0, surface = 0.0;
    Array_double_4 per_atom_surface;
    std::fill(per_atom_surface.begin(), per_atom_surface.end(), 0.0);

    // General idea: consider 24 subsimplexes with three orthogonal edges 
    // by constructing normals from weighted circumcenter to the faces and then to the edges
    for (int i = 0; i < 4; i++) {                             // iteration over cell facets
        const Point *p_vertex_i = pts[i];
        pts[i] = &V;     // replace vertex(i) in array  by weighted circumcenter of cell
        CGAL::Sign sV = orientation(*pts[0], *pts[1], *pts[2], *pts[3]);  // do V and vertex(i) on the same side of facet?

        const int a[3] = { (i+1)%4, (i+2)%4, (i+3)%4 };    // indicies of atoms of face opposite to cell vertex i
        const Point E = K::Plane_3(*pts[a[0]], *pts[a[1]], *pts[a[2]]).projection(V);  // projection of Voronoi vertex to facet
        const double z0 = sqrt((V-E).squared_length());    // length of normal to facet

        for (int j = 0; j < 3; j++) {                      // iteration over facet edges
            const int e[2] = { a[(j+1)%3], a[(j+2)%3] };   // indices of two edge atoms
            const Weighted_point WA[2] = { c->vertex(e[0])->point(), c->vertex(e[1])->point() };
            const Point A[2] = { WA[0].point(), WA[1].point() };
            const CGAL::Sign sE = coplanar_orientation(A[0], A[1], *pts[a[j]], E);

            const Point B = K::Line_3(A[0], A[1]).projection(E);  // projection of E to edge
            const double y0 = sqrt((E-B).squared_length());

            for (int k = 0; k < 2; k++) {                         // iteration over edge ends (two atoms)
                CGAL::Sign sB = CGAL::sign((B - A[k])*(A[(k+1)%2] - A[k]));  // negative if point B is in other direction than opposite A[(k+1)%2] looking from A[k]
                CGAL::Sign sF = sV * sE * sB;                                // final subsimplex sign
                double x0 = sqrt((B-A[k]).squared_length());
                double area;
                volume += sF * subsimplex_void_volume(x0, y0, z0, WA[k].weight(), &area);  // sum signed volume contribution from subsimplexes
                surface += sF * area;
                per_atom_surface[e[k]] += sF * area;
            }
        }

        pts[i] = p_vertex_i;  // place vertex(i) back to array
    }

    out_atom_surf = per_atom_surface;
    out_surf = surface;
    return volume;
}

// Returns whether the facet facet is attached to the cell
// that is used to represent facet
bool triangle_attached_to(const DT_3& dt, const Facet& facet) {
  Cell_handle cell = facet.first;
  int index = facet.second;

  if(dt.is_infinite(cell)) {
    return false;
  }

  Vertex_handle v1 = cell->vertex((index+1)%4);
  Vertex_handle v2 = cell->vertex((index+2)%4);
  Vertex_handle v3 = cell->vertex((index+3)%4);

  Vertex_handle w = cell->vertex(facet.second); 

  return CGAL::side_of_bounded_sphere(v1->point(),v2->point(),v3->point(),w->point())==CGAL::ON_BOUNDED_SIDE;
}

void draw_line_tetra(Delaunay &Tr,std::vector<Cell_handle> &cells,CGAL::Geomview_stream &gv) {
  if(! gv_on) return;
  CGAL::Color colors[] = {CGAL::BLUE,CGAL::GREEN,CGAL::YELLOW,CGAL::DEEPBLUE,
			  CGAL::PURPLE,CGAL::VIOLET,CGAL::ORANGE,CGAL::RED};
  int k = 0;
  for (std::vector<Cell_handle>::iterator it = cells.begin();it != cells.end();it++,k++) {
    Cell_handle c = *it;
    gv << colors[k % 8];
    Segment segs[4];
    for(int i = 0;i < 4;i++) {
      for(int j = i + 1;j < 4;j++) {
	Segment s = Segment(c->vertex(i)->point(),c->vertex(j)->point());
	gv << s;
      }
    }
  }
}

  void print(const Cell_handle c, const std::string& name) const
  {
    std::cerr << name << ":[" ;

    for ( int i=0; i<4 ; ++i )
      std::cerr << "[" << c->vertex(i)->point() << "]";

    std::cerr << "] ";
  }

Vertex_list fromCell(const Cell_handle& ch)
{
	Vertex_list the_list;
	for (auto i = 0; i < 4; i++)
	{
		the_list.push_back(ch->vertex(i));
	}
	return the_list;
}