Find boundary cells of C3t3 mesh

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

Find boundary cells of C3t3 mesh

Valentin Boltach
Good afternoon!

I ran into a little problem: Is there any interface for getting information
about the boundary cells of the figure's mesh?

CGAL has functions that can retrieve all information about boundary vertices
and faces but there's no info about cells.

I've tried one in some case incorrect variant, but it also doesn't work for
non-convex figure:

*0. We create a map from the vertices in the form of a key and the vertex
number as a value
*
    typedef typename Tr::Point Point_3;

    std::map<Point_3, int> V;
    Tr t = grid.triangulation();

    int k = 0;
    for (Tr::All_vertices_iterator it = t.all_vertices_begin(); it !=
t.all_vertices_end(); ++it) {
        if (k == 0) {
            k++;
            continue;
        }

        V[it->point()] = k;
        ++k;
    }

*1. We get a vector of all boundary faces (thanks to the method for output
mesh in .off format)
*
    C3t3::Subdomain_index sd_index = 0;
    bool normals_point_outside_of_the_subdomain = false;

    typedef std::vector<std::size_t> Face;

    std::vector<Point_3> points;
    std::vector<Face> faces;

    CGAL::Mesh_3::internal::facets_in_complex_3_to_triangle_soup(grid,
sd_index, points, faces,
                                                               
normals_point_outside_of_the_subdomain);


*2. We form a list of all possible surfaces for each cell and form a map,
the key of which are 3 vertices of the face, and the key is the vector of
cell numbers.
*
  std::map<std::set&lt;int>, vector<int>> cellNodes;
    typedef typename C3t3::Cell_iterator Cell_iterator;
    ofstream myfile;
    int k = 1;
    int facesPositions[4][3] = {{0,1,2}, {1,2,3}, {0,1,3}, {0,2,3}};
    for (Cell_iterator it = grid.cells_in_complex_begin(); it !=
grid.cells_in_complex_end(); ++it) {
        const Tr::Cell c(*it);

        for(int i = 0; i < 4; i++) {
            std::set<int> enodes;
            enodes.clear();
            for(int j = 0; j < 3; j++) {
                enodes.insert(V[c.vertex(facesPositions[i][j])->point()]);
            }
            cellNodes[enodes].push_back(k);
        }
        k++;
    }


*3. We iterate over all complex faces. If there is a match, then print the
number of the cell found
*
    size_t elemNumber = 0;

    typedef typename Tr::Cell_handle                                      
Cell_handle;
    typedef typename Tr::Vertex_handle                                    
Vertex_handle;
    typedef typename Tr::Weighted_point                                  
Weighted_point;
    typedef CGAL::Hash_handles_with_or_without_timestamps                
Hash_fct;

    typedef boost::unordered_map<Vertex_handle, std::size_t, Hash_fct>    
VHmap;
    typedef typename C3t3::Facets_in_complex_iterator                    
Ficit;

    VHmap vh_to_ids;
    std::size_t inum = 0;

    for(Ficit fit = grid.facets_in_complex_begin(),
                end = grid.facets_in_complex_end(); fit != end; ++fit)
    {
        Cell_handle c = fit->first;
        int s = fit->second;
        Face f;
        CGAL::Mesh_3::internal::resize(f, 3);

        typename C3t3::Subdomain_index cell_sdi = grid.subdomain_index(c);
        typename C3t3::Subdomain_index opp_sdi =
grid.subdomain_index(c->neighbor(s));

        if(cell_sdi != sd_index && opp_sdi != sd_index)
            continue;

        std::set<int> enodes;

        for(std::size_t i=1; i<4; ++i)
        {
            bool is_new;
            typename VHmap::iterator map_entry;
            Vertex_handle v = c->vertex((s+i)&3);
            enodes.insert(V[v->point()]);
            CGAL_assertion(v != Vertex_handle() &&
!grid.triangulation().is_infinite(v));
            boost::tie(map_entry, is_new) =
vh_to_ids.insert(std::make_pair(v, inum));
            if(is_new)
            {
                ++inum;
            }

            f[i-1] = map_entry->second;
        }

        if(((cell_sdi == sd_index) == (s%2 == 1)) ==
normals_point_outside_of_the_subdomain)
                std::swap(f[0], f[1]);

        const std::size_t fs = f.size();

        if(cellNodes.find(enodes) == cellNodes.end()) {
            LOG(INFO) << "NOT FOUND! " << enodes;
        }
        else {
            std::vector<int> elementsNumbers = cellNodes[enodes];
            LOG(INFO) << "Found! Cell#" << elementsNumbers.front();
        }
}

The problem is that Cell_handle contains information about *only 3 of its
vertices, not 4*, that's why I have to come up with options for iterating
over faces. How to be in this situation?




--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Find boundary cells of C3t3 mesh

Jane Tournois
Hello,

you can find all the info you need by using directly the triangulation
stored in the C3t3's data structure.

You can iterate on the C3t3 cells as you did with

`for (Cell_iterator it = grid.cells_in_complex_begin(); it
!=grid.cells_in_complex_end(); ++it)`

Then, for each cell, you can iterate on its 4 facets (see [1], Facet is
pair<Cell_handle, int>). For each facet, use `c3t3.is_in_complex(f)`
(see [2]) to check whether the facet is of interest for you.

Then, the function `mirror_facet(f)` [3] gives you the facet seen from
the adjacent cell, on the "other side" of the facet. So you have all
info about both cells adjacent to the complex facet.

Hope it helps!

Jane.

--
Jane Tournois, PhD
R&D Engineer at GeometryFactory
http://www.geometryfactory.com/

[1]
https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html#ad6a20b45e66dfb690bfcdb8438e9fcae

[2]
https://doc.cgal.org/latest/Mesh_3/classMeshComplex__3InTriangulation__3.html#a34e31b26c5578ab03a212f21e3ad703b

[3]
https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html#af386d55822ab1d47703f4500d642a7f2


Le 19/08/2019 à 22:50, Valentin Boltach a écrit :

> Good afternoon!
>
> I ran into a little problem: Is there any interface for getting information
> about the boundary cells of the figure's mesh?
>
> CGAL has functions that can retrieve all information about boundary vertices
> and faces but there's no info about cells.
>
> I've tried one in some case incorrect variant, but it also doesn't work for
> non-convex figure:
>
> *0. We create a map from the vertices in the form of a key and the vertex
> number as a value
> *
>      typedef typename Tr::Point Point_3;
>
>      std::map<Point_3, int> V;
>      Tr t = grid.triangulation();
>
>      int k = 0;
>      for (Tr::All_vertices_iterator it = t.all_vertices_begin(); it !=
> t.all_vertices_end(); ++it) {
>          if (k == 0) {
>              k++;
>              continue;
>          }
>
>          V[it->point()] = k;
>          ++k;
>      }
>
> *1. We get a vector of all boundary faces (thanks to the method for output
> mesh in .off format)
> *
>      C3t3::Subdomain_index sd_index = 0;
>      bool normals_point_outside_of_the_subdomain = false;
>
>      typedef std::vector<std::size_t> Face;
>
>      std::vector<Point_3> points;
>      std::vector<Face> faces;
>
>      CGAL::Mesh_3::internal::facets_in_complex_3_to_triangle_soup(grid,
> sd_index, points, faces,
>                                                                  
> normals_point_outside_of_the_subdomain);
>
>
> *2. We form a list of all possible surfaces for each cell and form a map,
> the key of which are 3 vertices of the face, and the key is the vector of
> cell numbers.
> *
>    std::map<std::set&lt;int>, vector<int>> cellNodes;
>      typedef typename C3t3::Cell_iterator Cell_iterator;
>      ofstream myfile;
>      int k = 1;
>      int facesPositions[4][3] = {{0,1,2}, {1,2,3}, {0,1,3}, {0,2,3}};
>      for (Cell_iterator it = grid.cells_in_complex_begin(); it !=
> grid.cells_in_complex_end(); ++it) {
>          const Tr::Cell c(*it);
>
>          for(int i = 0; i < 4; i++) {
>              std::set<int> enodes;
>              enodes.clear();
>              for(int j = 0; j < 3; j++) {
>                  enodes.insert(V[c.vertex(facesPositions[i][j])->point()]);
>              }
>              cellNodes[enodes].push_back(k);
>          }
>          k++;
>      }
>
>
> *3. We iterate over all complex faces. If there is a match, then print the
> number of the cell found
> *
>      size_t elemNumber = 0;
>
>      typedef typename Tr::Cell_handle
> Cell_handle;
>      typedef typename Tr::Vertex_handle
> Vertex_handle;
>      typedef typename Tr::Weighted_point
> Weighted_point;
>      typedef CGAL::Hash_handles_with_or_without_timestamps
> Hash_fct;
>
>      typedef boost::unordered_map<Vertex_handle, std::size_t, Hash_fct>
> VHmap;
>      typedef typename C3t3::Facets_in_complex_iterator
> Ficit;
>
>      VHmap vh_to_ids;
>      std::size_t inum = 0;
>
>      for(Ficit fit = grid.facets_in_complex_begin(),
>                  end = grid.facets_in_complex_end(); fit != end; ++fit)
>      {
>          Cell_handle c = fit->first;
>          int s = fit->second;
>          Face f;
>          CGAL::Mesh_3::internal::resize(f, 3);
>
>          typename C3t3::Subdomain_index cell_sdi = grid.subdomain_index(c);
>          typename C3t3::Subdomain_index opp_sdi =
> grid.subdomain_index(c->neighbor(s));
>
>          if(cell_sdi != sd_index && opp_sdi != sd_index)
>              continue;
>
>          std::set<int> enodes;
>
>          for(std::size_t i=1; i<4; ++i)
>          {
>              bool is_new;
>              typename VHmap::iterator map_entry;
>              Vertex_handle v = c->vertex((s+i)&3);
>              enodes.insert(V[v->point()]);
>              CGAL_assertion(v != Vertex_handle() &&
> !grid.triangulation().is_infinite(v));
>              boost::tie(map_entry, is_new) =
> vh_to_ids.insert(std::make_pair(v, inum));
>              if(is_new)
>              {
>                  ++inum;
>              }
>
>              f[i-1] = map_entry->second;
>          }
>
>          if(((cell_sdi == sd_index) == (s%2 == 1)) ==
> normals_point_outside_of_the_subdomain)
>                  std::swap(f[0], f[1]);
>
>          const std::size_t fs = f.size();
>
>          if(cellNodes.find(enodes) == cellNodes.end()) {
>              LOG(INFO) << "NOT FOUND! " << enodes;
>          }
>          else {
>              std::vector<int> elementsNumbers = cellNodes[enodes];
>              LOG(INFO) << "Found! Cell#" << elementsNumbers.front();
>          }
> }
>
> The problem is that Cell_handle contains information about *only 3 of its
> vertices, not 4*, that's why I have to come up with options for iterating
> over faces. How to be in this situation?
>
>
>
>
> --
> Sent from: http://cgal-discuss.949826.n4.nabble.com/
>

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Find boundary cells of C3t3 mesh

Valentin Boltach
Hi,
Thanks a lot for this answer. I understood everything except iteration on
facets.
Because i'm new in cgal, it's quite difficult to understand all interfaces
:)
Can you describe this moment in detail?



--
Sent from: http://cgal-discuss.949826.n4.nabble.com/

--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss


Reply | Threaded
Open this post in threaded view
|

Re: Find boundary cells of C3t3 mesh

Jane Tournois
Hi,

for each `Cell_handle c`, and i an integer between 0 and 3, Facet(c, i)
is the ith facet of c, seen from inside c. So you only need to iterate on i.

All the data structure is explained in the
Triangulation_data_structure_3 documentation :

https://doc.cgal.org/latest/TDS_3/classTriangulationDataStructure__3.html


Le 23/08/2019 à 22:54, Valentin Boltach a écrit :

> Hi,
> Thanks a lot for this answer. I understood everything except iteration on
> facets.
> Because i'm new in cgal, it's quite difficult to understand all interfaces
> :)
> Can you describe this moment in detail?
>
>
>
> --
> Sent from: http://cgal-discuss.949826.n4.nabble.com/
>
--
Jane Tournois, PhD
R&D Engineer at GeometryFactory
http://www.geometryfactory.com/


--
You are currently subscribed to cgal-discuss.
To unsubscribe or access the archives, go to
https://sympa.inria.fr/sympa/info/cgal-discuss