C++ Triangulation::finite_vertices_begin方法代码示例

void Foam::DelaunayMeshTools::writeFixedPoints
(
    const fileName& fName,
    const Triangulation& t
)
{
    OFstream str(fName);

    Pout<< nl
        << "Writing fixed points to " << str.name() << endl;

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            t.finite_vertices_begin();
        vit != t.finite_vertices_end();
        ++vit
    )
    {
        if (vit->fixed())
        {
            meshTools::writeOBJ(str, topoint(vit->point()));
        }
    }
}

void solve(Triangulation DT, Segment rect[],vector<Point> points)
{
	map<TP,double> vertex_dist;
	for (vertex_iterator vit = DT.finite_vertices_begin();vit!=DT.finite_vertices_end();vit++)
	{
		TP vertex_p = vit->point();
		double min_dist = min_segment(rect, vertex_p);
		vertex_dist[vertex_p] = min_dist;
	}	
	
	for(Edge_iterator eit=	DT.finite_edges_begin(); eit!=DT.finite_edges_end();eit++)
	{
		double min_dist = CGAL::to_double(DT.segment(eit).squared_length()/4);
		TP vertex_p1 = eit->first->vertex(DT.cw(eit->second))->point();
		TP vertex_p2 = eit->first->vertex(DT.ccw(eit->second))->point();

		vertex_dist[vertex_p1] = std::min(vertex_dist[vertex_p1], min_dist);
		vertex_dist[vertex_p2] = std::min(vertex_dist[vertex_p2], min_dist);
	}

	vector<double> dist;
	for (std::map<TP,double>::iterator it = vertex_dist.begin();it!= vertex_dist.end(); it++)
	{
		dist.push_back(it->second);
	}

	sort(dist.begin(),dist.end());

	int f,m,l;
	f = time(dist[0]);
	m = time(dist[points.size()/2]);
	l = time(dist[points.size() - 1]);

	cout<<f<<" "<<m<<" "<<l<<endl;	
}

Foam::tmp<Foam::Field<Type>> filterFarPoints
(
    const Triangulation& mesh,
    const Field<Type>& field
)
{
    tmp<Field<Type>> tNewField(new Field<Type>(field.size()));
    Field<Type>& newField = tNewField.ref();

    label added = 0;
    label count = 0;

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            mesh.finite_vertices_begin();
        vit != mesh.finite_vertices_end();
        ++vit
    )
    {
        if (vit->real())
        {
            newField[added++] = field[count];
        }

        count++;
    }

    newField.resize(added);

    return tNewField;
}

void u_new(Triangulation& T, const FT dt ) {

  for(F_v_it fv=T.finite_vertices_begin();
      fv!=T.finite_vertices_end();
      fv++) {

//  for(F_v_it fv=Tp.finite_vertices_begin();
//      fv!=Tp.finite_vertices_end();
//      fv++) {

    Vector_2 Ustar = fv->Ustar.val() ;
    Vector_2 gradp = fv->gradp.val() ;
    Vector_2 U = Ustar - dt * gradp;

    // relaxation mixing .-
    FT alpha=simu.alpha();
    Vector_2 U0=fv->U() ;
    Vector_2 U_mix = alpha*U0+ (1-alpha)*U ;

    fv->U.set( U_mix );
    fv->Delta_U.set(  U_mix - fv->Uold.val()  );
    
  }

  return;

}

void update_half_velocity( Triangulation& Tp , const bool overdamped ) {

  if (overdamped) return;
  
   for(F_v_it fv=Tp.finite_vertices_begin();
       fv!=Tp.finite_vertices_end();
       fv++) {

    Vector_2  v  = fv->U();

    // if (overdamped) 
    //    fv->U.set( v );
    // else {

    Vector_2  v0 = fv->Uold();
      //    Vector_2  v_star = fv->Ustar();

    fv->U.set(  2 * v - v0 );

      //   fv->U.set(  v + v_star - v0 );

      //    }

  }
  
  return;

}

void Foam::DelaunayMeshTools::writeInternalDelaunayVertices
(
    const fileName& instance,
    const Triangulation& t
)
{
    pointField internalDelaunayVertices(t.number_of_vertices());

    label vertI = 0;

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            t.finite_vertices_begin();
        vit != t.finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalPoint())
        {
            internalDelaunayVertices[vertI++] = topoint(vit->point());
        }
    }

    internalDelaunayVertices.setSize(vertI);

    pointIOField internalDVs
    (
        IOobject
        (
            "internalDelaunayVertices",
            instance,
            t.time(),
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        internalDelaunayVertices
    );

    Info<< nl
        << "Writing " << internalDVs.name()
        << " to " << internalDVs.instance()
        << endl;

    internalDVs.write();
}

void update_half_alpha( Triangulation& Tp ) {

   for(F_v_it fv=Tp.finite_vertices_begin();
       fv!=Tp.finite_vertices_end();
       fv++) {

    FT al  = fv->alpha();

    FT al0 = fv->alpha0();

    fv->alpha.set(  2 * al - al0 );

  }

  return;

}

// ---------------------------------------------------------
// initialize
// ----------
// Initialize some of the attributes of the triangulation.
// ---------------------------------------------------------
void
initialize(Triangulation &triang)
{
   // set vertex id.
   int id = 0;
   for(FVI vit = triang.finite_vertices_begin();
	vit != triang.finite_vertices_end(); vit ++)
   {
	   vit->id = id++;
	   vit->visited = false;
    	   vit->bad = false;
    	   vit->bad_neighbor = false;
   }

   // set cell id.
   id = 0;
   for(ACI cit = triang.all_cells_begin();
	cit != triang.all_cells_end(); cit ++)
   {
	   cit->id = id++;
    	   cit->visited = false;
    	   cit->outside = false;
    	   cit->transp = false;

	   for(int id = 0 ; id < 4; id++)
	   {
      		cit->set_cocone_flag(id,false);
      		cit->neighbor(id)->set_cocone_flag(cit->neighbor(id)->index(cit),false);
      		cit->bdy[id] = false;
      		cit->opaque[id] = false;
      		for(int k = 0; k < 4; k ++)
	    		cit->umbrella_member[id][k] = -1;
    	   }
    	   
	   // set the convex hull points.
    	   if(! triang.is_infinite(cit)) continue;

    	   for(int i = 0; i < 4; i ++)
    	   {
	    	   if(! triang.is_infinite(cit->vertex(i))) continue;
	    	   cit->vertex((i+1)%4)->set_convex_hull(true);
	    	   cit->vertex((i+2)%4)->set_convex_hull(true);
	    	   cit->vertex((i+3)%4)->set_convex_hull(true);
    	   }
   }
}

void move_info(Triangulation& T) {

  for(F_v_it fv=T.finite_vertices_begin();
      fv!=T.finite_vertices_end();
      fv++) {

    Periodic_point rr=T.periodic_point(fv);

    Point r=T.point(rr); // current point

    fv->alpha0.set( fv->alpha() );
    fv->rold.set(r);
    fv->Uold.set(fv->U());
    //    fv->Ustar.set(fv->U());
    //    fv->p.set(fv->p() +  fv->pstar() );
  }


  return;
}

Foam::tmp<Foam::pointField> Foam::DelaunayMeshTools::allPoints
(
    const Triangulation& t
)
{
    tmp<pointField> tpts(new pointField(t.vertexCount(), point::max));
    pointField& pts = tpts.ref();

    for
    (
        typename Triangulation::Finite_vertices_iterator vit =
            t.finite_vertices_begin();
        vit != t.finite_vertices_end();
        ++vit
    )
    {
        if (vit->internalOrBoundaryPoint() && !vit->referred())
        {
            pts[vit->index()] = topoint(vit->point());
        }
    }

    return tpts;
}

void u_star(Triangulation& T, FT dt , bool semi ) { 

  for(F_v_it fv=T.finite_vertices_begin();
      fv!=T.finite_vertices_end();
      fv++) {

    Vector_2 U0 = fv->Uold() ;

    Vector_2 f;

    if(semi)
      f = simu.mu() * fv->laplU() - fv->gradp() + fv->force();
    else
      f = simu.mu() * fv->laplU() + fv->force(); // Include all forces here, but for the pressure grad!

    Vector_2 Ustar = U0 + dt * f;

    fv->Ustar.set(Ustar) ;

  }

  return;

}

FT move(Triangulation& Tp, const FT dt , FT& dd0 ) {

  vector<data_kept> prev;

  FT dd2=0;

  bool first=false;  // debug

  for(F_v_it fv=Tp.finite_vertices_begin();
      fv!=Tp.finite_vertices_end();
      fv++) {
    data_kept data(fv);

    Vector_2  vel = fv->U();

    Vector_2 disp = dt * vel;

    Periodic_point rr=Tp.periodic_point(fv);

    Point rnow=Tp.point(rr); // current point

    Point r0=fv->rold(); // starting point

    Point rnew= r0 + disp;

    Vector_2 disp2 = per_vect(rnew,rnow);

    FT rel_disp = sqrt(disp2.squared_length() ) / simu.h();

    FT rel_disp0= sqrt( disp.squared_length() ) / simu.h();

    if(first) {
     cout
       << "r0 " << r0 << "  "
       << "rnow " << rnow << "  "
       << "rnew " << rnew << "  "
       << "disp2 " << disp2 << "  "
       << " idx " << fv->idx() << " "
       << "rel_disp " << rel_disp
       << endl ;
     first=false;
    }

    dd2 += rel_disp;

    dd0 += rel_disp0;

    //    cout << "New position: " << r0 ;

    data.pos = per_point( rnew );

    //    cout << " ---> " << data.pos  << endl ;

    prev.push_back (data);

  }

//  cout << "relative displacement " << sqrt(dd2)/simu.no_of_points()/simu.h()  << endl ;
  dd2 /= simu.no_of_particles();
  dd0 /= simu.no_of_particles();

  //  cout << "relative displacement " << dd2 << endl ;

  Tp.clear(); // clears the triangulation !!

  for(vector<data_kept>::iterator data=prev.begin();
      data!=prev.end();
      data++) {

    //    cout << "Inserting back at " << data->pos << endl ;

    Vertex_handle fv=Tp.insert(data->pos);

    data->restore(fv);

    // return info to vertices


  }

//  cout << "Insertion done" << endl ;
  Tp.convert_to_1_sheeted_covering();

  return dd2;
}

int main() {

    while(true) {
        int bacteria_count;
        cin >> bacteria_count;

        // kill switch for application
        if(bacteria_count == 0) {
            break;
        }

        // read in boundaries of the dish
        double left_border, right_border, bottom_border, top_border;
        cin >> left_border >> bottom_border >> right_border >> top_border;

        // collect bacteria's center information
        vector<K::Point_2> bacteria_centers;
        bacteria_centers.reserve(bacteria_count);
        for(int i = 0; i < bacteria_count; i++) {
            double bacteria_x, bacteria_y;
            cin >> bacteria_x >> bacteria_y;

            bacteria_centers.push_back(K::Point_2(bacteria_x, bacteria_y));
        }

        // create triangulation
        Triangulation triang;
        triang.insert(bacteria_centers.begin(), bacteria_centers.end());

        // keep track of the distances for each bacteria
        map<Triangulation::Point, double> distances;
        //distances.reserve(bacteria_count);

        // calculate initial distance: distance between the bacteria and the nearest dish boundary
        for(Triangulation::Finite_vertices_iterator vertex_iter = triang.finite_vertices_begin(); vertex_iter != triang.finite_vertices_end(); ++vertex_iter) {
            Triangulation::Point vertex = vertex_iter->point();
            distances[vertex] = min(
                min(vertex.x() - left_border, right_border - vertex.x()), // left/right minimum
                min(vertex.y() - bottom_border, top_border - vertex.y()) // top/bottom minimum
            );

            distances[vertex] *= distances[vertex]; // square distance as we work with squared ones
        }

        // compute distance to other two neighbours and update distance if it is smaller
        for(Triangulation::Finite_edges_iterator edge_iter = triang.finite_edges_begin(); edge_iter != triang.finite_edges_end(); ++edge_iter) {
            Triangulation::Vertex_handle vertex1 = edge_iter->first->vertex(triang.cw(edge_iter->second));
            Triangulation::Vertex_handle vertex2 = edge_iter->first->vertex(triang.ccw(edge_iter->second));

            Triangulation::Point vertex1_point = vertex1->point();
            Triangulation::Point vertex2_point = vertex2->point();

            // calculate distance of the points of both vertex and half them (divide by 4 as distance is squared and 4 = 2^2)
            double vertex_distance = CGAL::to_double(CGAL::squared_distance(vertex1_point, vertex2_point)) / 4;

            // update distances to minimum
            distances[vertex1_point] = min(distances[vertex1_point], vertex_distance);
            distances[vertex2_point] = min(distances[vertex2_point], vertex_distance);
        }

        // now we know the minimum distance for each bacteria to another one or the borders of the dish

        // extract distances into a vector and sort it
        vector<double> only_distances;
        only_distances.reserve(bacteria_count);
        for(map<Triangulation::Point, double>::iterator iter = distances.begin(); iter != distances.end(); ++iter) {
            only_distances.push_back(iter->second);
        }

        // sort distances
        sort(only_distances.begin(), only_distances.end());

        // print out information
        cout << hours(only_distances[0]) << " " << hours(only_distances[bacteria_count/2]) << " " << hours(only_distances[bacteria_count - 1]) << endl;
    }
}

bool Foam::conformalVoronoiMesh::distributeBackground(const Triangulation& mesh)
{
    if (!Pstream::parRun())
    {
        return false;
    }

    Info<< nl << "Redistributing points" << endl;

    timeCheck("Before distribute");

    label iteration = 0;

    scalar previousLoadUnbalance = 0;

    while (true)
    {
        scalar maxLoadUnbalance = mesh.calculateLoadUnbalance();

        if
        (
            maxLoadUnbalance <= foamyHexMeshControls().maxLoadUnbalance()
         || maxLoadUnbalance <= previousLoadUnbalance
        )
        {
            // If this is the first iteration, return false, if it was a
            // subsequent one, return true;
            return iteration != 0;
        }

        previousLoadUnbalance = maxLoadUnbalance;

        Info<< "    Total number of vertices before redistribution "
            << returnReduce(label(mesh.number_of_vertices()), sumOp<label>())
            << endl;

        const fvMesh& bMesh = decomposition_().mesh();

        volScalarField cellWeights
        (
            IOobject
            (
                "cellWeights",
                bMesh.time().timeName(),
                bMesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            bMesh,
            dimensionedScalar("weight", dimless, 1e-2),
            zeroGradientFvPatchScalarField::typeName
        );

        meshSearch cellSearch(bMesh, polyMesh::FACE_PLANES);

        labelList cellVertices(bMesh.nCells(), label(0));

        for
        (
            typename Triangulation::Finite_vertices_iterator vit
                = mesh.finite_vertices_begin();
            vit != mesh.finite_vertices_end();
            ++vit
        )
        {
            // Only store real vertices that are not feature vertices
            if (vit->real() && !vit->featurePoint())
            {
                pointFromPoint v = topoint(vit->point());

                label cellI = cellSearch.findCell(v);

                if (cellI == -1)
                {
//                     Pout<< "findCell conformalVoronoiMesh::distribute "
//                         << "findCell "
//                         << vit->type() << " "
//                         << vit->index() << " "
//                         << v << " "
//                         << cellI
//                         << " find nearest cellI ";

                    cellI = cellSearch.findNearestCell(v);
                }

                cellVertices[cellI]++;
            }
        }

        forAll(cellVertices, cI)
        {
            // Give a small but finite weight for empty cells.  Some
            // decomposition methods have difficulty with integer overflows in
            // the sum of the normalised weight field.
            cellWeights.internalField()[cI] = max
            (
                cellVertices[cI],
                1e-2
            );
        }
//.........这里部分代码省略.........