multiple copy of neighbors in periodic 2d triangulation

classic Classic list List threaded Threaded
2 messages Options
Reply | Threaded
Open this post in threaded view
|

multiple copy of neighbors in periodic 2d triangulation

Roya Rahmati
Hello,
I am using periodic Delaunay triangulation in my code, and producing for each vertex all neighboring vertices. For this I use Edge iterator, since in my case it will be much more faster than Vertex iterator.

But, in some of my simulations, especially for particle numbers<1000 instead of one special neighbor for each vertex I get 8 or 9 or ... copy of the same neighbor. For example I get something like this for the vertex 0:  Nei[0]=[2,2,2,2,2,2,2,2,4,4,4,...]

I couldn't find an example using edge iterator in CGAL website, and nothing related in stackoverflow. I am wondering whether this implementation is correct? What should I add to my code in order to get one copy per each neighbor, I can use STL::set, but I would like to know the source of problem.
A possible work around would be to replicate main simulation box, and use non-periodic version.
here is the code snippet:
    typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
    typedef CGAL::Periodic_2_triangulation_traits_2<Kernel> Gt;
    typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, Gt> Vb;
    typedef CGAL::Periodic_2_triangulation_face_base_2<Gt> Fb;
    typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
    typedef CGAL::Periodic_2_Delaunay_triangulation_2<Gt, Tds> Triangulation;
    typedef Triangulation::Iso_rectangle Iso_rectangle;
    typedef Triangulation::Edge_iterator Edge_iterator;
    typedef Triangulation::Vertex_handle Vertex_handle;
    typedef Triangulation::Point Point;
    typedef vector<pair<Point, unsigned> > Vector_Paired;
    Vector_Paired points;
    Iso_rectangle domain(0,0,L,L);

    for(int iat = 0; iat < N; iat++)
        {
     points.push_back(make_pair(Point(r_tot[iat][0],r_tot[iat][1]),iat));
        }
    Triangulation T(points.begin(), points.end(), domain);

    for(Edge_iterator ei=T.finite_edges_begin(); ei!=T.finite_edges_end(); ei++)
        {
         
          Triangulation::Face& f = *(ei->first);
          int ii = ei->second;
          Vertex_handle vi = f.vertex(f.cw(ii));     
          Vertex_handle vj = f.vertex(f.ccw(ii));    
          int iat = vi->info();
          int jat = vj->info();

          VecInd[iat].push_back(jat);
          VecInd[jat].push_back(iat);
             
        }
Regards,
Reply | Threaded
Open this post in threaded view
|

Re: multiple copy of neighbors in periodic 2d triangulation

MaelRL

If your point set is not geometrically well spaced, it's possible that the triangulation of these points do not form a simplicial complex over the flat torus (in other words, there are short cycles in the triangulation). In this case, the algorithm uses 8 copies of the triangulation to artificially create a simplicial complex. You can check if this is the case using the function "is_triangulation_in_1_sheet()" (https://github.com/CGAL/cgal/blob/master/Periodic_2_triangulation_2/include/CGAL/Periodic_2_triangulation_2.h#L4108) and read more about these mechanisms in the User Manual (https://doc.cgal.org/latest/Periodic_2_triangulation_2/index.html#title1).

When copies are being used, iterating over the edges will indeed give you exactly what the underlying data structure has : 9 entities for each edge. To get unique ones, you can simply filter 8 out of the 9 by looking at the offset of the vertices of the edge. This is what is done in the iterator that returns unique periodic segments. Unfortunately, you want edges and this iterator converts directly to the geometry of the edge (the segment). Nevertheless, you can simply use the main filtering function from that iterator, that is: "is_canonical()" (https://github.com/CGAL/cgal/blob/c68cf8fc4c850f8cd84c6900faa781286a7117ed/Periodic_2_triangulation_2/include/CGAL/Periodic_2_triangulation_iterators_2.h#L505). This function will look at the offset of the two vertices of your edges, and keep only those that have at least one vertex in the first copy of the domain, which is enough to make it unique.

Best,
Mael

On 2019-05-02 18:12, Roya Rahmati wrote:
Hello,
I am using periodic Delaunay triangulation in my code, and producing for each vertex all neighboring vertices. For this I use Edge iterator, since in my case it will be much more faster than Vertex iterator.

But, in some of my simulations, especially for particle numbers<1000 instead of one special neighbor for each vertex I get 8 or 9 or ... copy of the same neighbor. For example I get something like this for the vertex 0:  Nei[0]=[2,2,2,2,2,2,2,2,4,4,4,...]

I couldn't find an example using edge iterator in CGAL website, and nothing related in stackoverflow. I am wondering whether this implementation is correct? What should I add to my code in order to get one copy per each neighbor, I can use STL::set, but I would like to know the source of problem.
A possible work around would be to replicate main simulation box, and use non-periodic version.
here is the code snippet:
    typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
    typedef CGAL::Periodic_2_triangulation_traits_2<Kernel> Gt;
    typedef CGAL::Triangulation_vertex_base_with_info_2<unsigned int, Gt> Vb;
    typedef CGAL::Periodic_2_triangulation_face_base_2<Gt> Fb;
    typedef CGAL::Triangulation_data_structure_2<Vb, Fb> Tds;
    typedef CGAL::Periodic_2_Delaunay_triangulation_2<Gt, Tds> Triangulation;
    typedef Triangulation::Iso_rectangle Iso_rectangle;
    typedef Triangulation::Edge_iterator Edge_iterator;
    typedef Triangulation::Vertex_handle Vertex_handle;
    typedef Triangulation::Point Point;
    typedef vector<pair<Point, unsigned> > Vector_Paired;
    Vector_Paired points;
    Iso_rectangle domain(0,0,L,L);

    for(int iat = 0; iat < N; iat++)
        {
     points.push_back(make_pair(Point(r_tot[iat][0],r_tot[iat][1]),iat));
        }
    Triangulation T(points.begin(), points.end(), domain);

    for(Edge_iterator ei=T.finite_edges_begin(); ei!=T.finite_edges_end(); ei++)
        {
         
          Triangulation::Face& f = *(ei->first);
          int ii = ei->second;
          Vertex_handle vi = f.vertex(f.cw(ii));     
          Vertex_handle vj = f.vertex(f.ccw(ii));    
          int iat = vi->info();
          int jat = vj->info();

          VecInd[iat].push_back(jat);
          VecInd[jat].push_back(iat);
             
        }
Regards,