CGAL users discussion list

[<>]

Dear all,

I am trying to expand a periodic mesh and have it returned in a new c3t3 structure instead of having to save a .mesh file and load it afterwards. For this I have tried modifying the File_medit.h found in the Periodic_3_mesh_3/IO folder but am having a few issues. First of all, the tetrahedrals seem to not be properly added when I call the build_triangulation function although the finite_cell vector has all collected values in it. I am having a similar issue with the border_facets, although in this case I am also not sure what data I am supposed to place in surface_patch_id and if perhaps this is causing the issues. I have included the code I am using bellow. Any help is appreciated.

Regards, Kim

/* ==========================================================================

 *  Based on CGAL's File_medit.h. See CGAL Periodic_3_mesh_3/IO.

 * ==========================================================================*/

#ifndef CGAL_PERIODIC_REPLICATOR_H

#define CGAL_PERIODIC_REPLICATOR_H

#include <CGAL/license/Periodic_3_mesh_3.h>

#include <CGAL/Mesh_3/io_signature.h>

#include <CGAL/Mesh_3/tet_soup_to_c3t3.h>

//

#include <CGAL/array.h>

#include <CGAL/assertions.h>

#include <CGAL/IO/File_medit.h>

#include <boost/unordered_map.hpp>

#include <algorithm>

#include <fstream>

#include <iostream>

#include <cstring>

#include <sstream>

#include <cstring>

#include <map>

namespace CGAL {

  namespace periodic_replicator {

    /**

     * \brief expands a periodic mesh to the c3t3 variable

     * \param c3t3_periodic the mesh

     * \param occurrence_count the number of copies that are printed

     * \param distinguish_copies if set to `true`, each copy is assigned a unique color.

     *                           Otherwise, all domains are drawn with subdomain index-based colors.

     * \param rebind if set to `true`, labels of cells are rebinded into [1..nb_of_labels]

     * \param c3t3 the output mesh

     */

    template <class C3T3_Periodic, class C3T3, bool rebind, bool no_patch>  ////template <class C3T3_Periodic>

    void replicate_periodic_mesh(const C3T3_Periodic& c3t3_periodic,

                                int occurrence_count[3],

                                const bool distinguish_copies,

                                C3T3& c3t3)

    {

      #ifdef CGAL_MESH_3_IO_VERBOSE

        std::cerr << "Output to medit:\n";

      #endif

      CGAL_precondition(c3t3_periodic.triangulation().is_1_cover());

      typedef CGAL::Mesh_3::Medit_pmap_generator<C3T3_Periodic, rebind, no_patch>  Generator;

      typedef typename Generator::Cell_pmap                               Cell_pmap;

      typedef typename Generator::Facet_pmap                              Facet_pmap;

      typedef typename Generator::Facet_pmap_twice                        Facet_pmap_twice;

      typedef typename Generator::Vertex_pmap                             Vertex_pmap;

      //Generator().print_twice();

      Cell_pmap cell_pmap(c3t3_periodic);

      Facet_pmap facet_pmap(c3t3_periodic, cell_pmap);

      Facet_pmap_twice facet_pmap_twice(c3t3_periodic, cell_pmap);

      Vertex_pmap vertex_pmap(c3t3_periodic, cell_pmap, facet_pmap);

  //    const Facet_index_property_map_twice& facet_twice_pmap = Facet_index_property_map_twice(),

      const bool print_each_facet_twice = false;//Generator().print_twice();

     

      CGAL_precondition(occurrence_count[0] >= 1 && occurrence_count[1] >= 1 && occurrence_count[2] >= 1);

      // Periodic typedef

      typedef typename C3T3_Periodic::Triangulation  Triangulation_periodic;

      typedef Triangulation_periodic                 Tr_periodic;

      typedef typename Tr_periodic::Bare_point       Bare_point;

      typedef typename Tr_periodic::Weighted_point   Weighted_point;

      typedef typename C3T3_Periodic::Vertex_handle  Vertex_handle;

      typedef typename C3T3_Periodic::Cell_handle    Cell_handle;

      typedef typename Tr_periodic::Vertex_iterator  Vertex_iterator;

      typedef typename C3T3_Periodic::Facet_iterator Facet_iterator;

      typedef typename C3T3_Periodic::Cell_iterator  Cell_iterator;

      typedef typename Tr_periodic::Offset           Offset;

      const Triangulation_periodic& tr_periodic = c3t3_periodic.triangulation();

      // Non-periodic typedef

      typedef typename C3T3::Triangulation  Triangulation;

      typedef Triangulation                 Tr;

      typedef typename Tr::Point Point_3;

      typedef std::array<int, 3> Facet; // 3 = id

      typedef std::array<int, 5> Tet_with_ref; // first 4 = id, fifth = reference

      //

      std::vector<Point_3> points;

      //std::map<Facet, typename Tr::Cell::Surface_patch_index> border_facets_;

      std::map<Facet, typename C3T3::Surface_patch_index> border_facets;

      std::vector<Tet_with_ref> finite_cells;

      //

      int number_of_vertices = static_cast<int>(tr_periodic.number_of_vertices());

      int number_of_facets = static_cast<int>(c3t3_periodic.number_of_facets());

      int number_of_cells = static_cast<int>(c3t3_periodic.number_of_cells());

      // number of reproductions over the different axes

      int Ox_rn = ceil(occurrence_count[0]);

      int Oy_rn = ceil(occurrence_count[1]);

      int Oz_rn = ceil(occurrence_count[2]);

      int occ_mult = Ox_rn * Oy_rn * Oz_rn;

      #ifdef CGAL_PERIODIC_3_MESH_3_VERBOSE

        std::cerr << "Expanding periodic mesh... " << std::endl;

        std::cerr << "occurrences over each axis: "

                  << Ox_rn << " " << Oy_rn << " " << Oz_rn << std::endl;

        std::cerr << number_of_vertices << " vertices" << std::endl;

        std::cerr << number_of_facets << " facets" << std::endl;

        std::cerr << number_of_cells << " cells" << std::endl;

      #endif

      std::ofstream os("output_control.mesh");

      os << std::setprecision(17);

      // Vertices

      int medit_number_of_vertices = (Ox_rn + 1) * (Oy_rn + 1) * (Oz_rn + 1) * number_of_vertices;

      os << "MeshVersionFormatted 1\n"

         << "Dimension 3\n"

         << "Vertices\n" << medit_number_of_vertices << std::endl;

      // Build the set of points that is needed to draw all the elements.

      // On each axis, we repeat n+1 times the point, where 'n' is the number of

      // instances of the mesh that will be printed over that axis. This is because

      // a cell 'c' might have point(c,i) that is equal to v with an offset 2

      boost::unordered_map<Vertex_handle, int> V;

      int inum = 1; // '1' because medit ids start at 1

      for(int i=0; i<=Oz_rn; ++i) {

        for(int j=0; j<=Oy_rn; ++j) {

          for(int k=0; k<=Ox_rn; ++k) {

            for(Vertex_iterator vit = tr_periodic.vertices_begin(); vit != tr_periodic.vertices_end(); ++vit) {

              if(i == 0 && j == 0 && k == 0)

                V[vit] = inum++;

              const Offset off(k, j, i);

              const Weighted_point& p = tr_periodic.point(vit);

              const Bare_point bp = tr_periodic.construct_point(p, off);

              int id;

              if(i >= 1 || j >= 1 || k >= 1)

                id = 7;

              else

                id = tr_periodic.off_to_int(off);

              os << CGAL::to_double(bp.x()) << ' '

                 << CGAL::to_double(bp.y()) << ' '

                 << CGAL::to_double(bp.z()) << ' ';

              if(!distinguish_copies || occ_mult == 1)

                os << get(vertex_pmap, vit) << '\n';

              else

                os << 32 * id + 1 << '\n';

              points.push_back(Point_3(CGAL::to_double(bp.x()),

                                       CGAL::to_double(bp.y()),

                                       CGAL::to_double(bp.z())));

            }

          }

        }

      }

      // Triangles

      int medit_number_of_triangles = occ_mult * number_of_facets;

      if(print_each_facet_twice)

        medit_number_of_triangles *= 2;

      os << "Triangles\n" << medit_number_of_triangles << std::endl;

      for(int i=0; i<Oz_rn; ++i) {

        for(int j=0; j<Oy_rn; ++j) {

          for(int k=0; k<Ox_rn; ++k) {

            const Offset off(k, j, i);

            for(Facet_iterator fit = c3t3_periodic.facets_begin(); fit != c3t3_periodic.facets_end(); ++fit) {

              typename Tr::Cell::Surface_patch_index surface_patch_id;

              Facet facet;

              for(int l=0; l<4; ++l) {

                if(l == fit->second)

                  continue;

                Cell_handle c = fit->first;

                Vertex_handle v = c->vertex(l);

                const Offset combined_off = tr_periodic.combine_offsets(

                                              off, tr_periodic.int_to_off(c->offset(l)));

                const int vector_offset = combined_off.x() +

                                          combined_off.y() * (Ox_rn + 1) +

                                          combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);

                const int id = vector_offset * number_of_vertices + V[v];

                CGAL_assertion(1 <= id && id <= medit_number_of_vertices);

                os << id << " ";

                //facet[l] = id;

                int temp_ = id;

                facet[l] = temp_;

              }

              // For multiple copies, color to distinguish copies rather than to distinguish subdomains

              //int temptest;

              if(!distinguish_copies || occ_mult == 1){

                os << get(facet_pmap, *fit) << '\n';

              } else {

                os << 1 + k + 3*j + 9*i << '\n';

              }

              surface_patch_id.first = 0;

              surface_patch_id.second = 1;

              border_facets.insert(std::make_pair(facet, surface_patch_id));

              // Print triangle again if needed

              if(print_each_facet_twice) {

                for(int l=0; l<4; ++l) {

                  if(l == fit->second)

                    continue;

                  Cell_handle c = fit->first;

                  Vertex_handle v = c->vertex(l);

                  const Offset combined_off = tr_periodic.combine_offsets(

                                                off, tr_periodic.int_to_off(c->offset(l)));

                  const int vector_offset = combined_off.x() +

                                            combined_off.y() * (Ox_rn + 1) +

                                            combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);

                  const int id = vector_offset * number_of_vertices + V[v];

                  CGAL_assertion(1 <= id && id <= medit_number_of_vertices);

                  os << id << " ";

                  int temp_ = id;

                  facet[l] = temp_;

                }

                if(!distinguish_copies || occ_mult == 1)

                  os << get(facet_pmap_twice, *fit) << '\n';

                else

                  os << 1 + k + 3*j + 9*i << '\n';

                  //surface_patch_id = ???;

                  //border_facets.insert(std::make_pair(facet, surface_patch_id));

              }

            }

          }

        }

      }

      // Tetrahedra

      os << "Tetrahedra\n" << occ_mult * number_of_cells << std::endl;

      for(int i=0; i<Oz_rn; ++i) {

        for(int j=0; j<Oy_rn; ++j) {

          for(int k=0; k<Ox_rn; ++k) {

            const Offset off(k, j, i);

            for(Cell_iterator cit = c3t3_periodic.cells_begin(); cit !=c3t3_periodic.cells_end(); ++cit) {

              Tet_with_ref t;

              for(int l=0; l<4; ++l) {

                const Offset combined_off = tr_periodic.combine_offsets(

                                              off, tr_periodic.int_to_off(cit->offset(l)));

                const int vector_offset = combined_off.x() +

                                          combined_off.y() * (Ox_rn + 1) +

                                          combined_off.z() * (Ox_rn + 1) * (Oy_rn + 1);

                const int id = vector_offset * number_of_vertices + V[cit->vertex(l)];

                CGAL_assertion(1 <= id && id <= medit_number_of_vertices);

                os << id << " ";

                t[l] = id - 1;

              }

                   

              // For multiple copies, color to distinguish copies rather than to distinguish subdomains

              if(!distinguish_copies || occ_mult == 1) {

                os << get(cell_pmap, cit) << '\n';

                t[4] = get(cell_pmap, cit);

              } else {

                os << 1 + k + 3*j + 9*i << '\n';

                t[4] = 1 + k + 3*j + 9*i;

              }

              finite_cells.push_back(t);

            }

          }

        }

      }

      os << "End" << std::endl;

      //

      //std::vector<typename Tr::Vertex_handle> vertices(0);

      std::vector<typename Tr::Vertex_handle> vertices(points.size() + 1);

      bool is_well_built = build_triangulation<Tr, true>(c3t3.triangulation(), points, finite_cells, border_facets, vertices, true/*verbose*/, false/*replace_domain*/);

      // Check

      std::cout << points.size() << " points" << std::endl;

      std::cout << border_facets.size() << " border facets" << std::endl;

      std::cout << finite_cells.size() << " cells" << std::endl;

      std::cout << c3t3.number_of_facets() << " border facets expanded c3t3" << std::endl;

      std::cout << c3t3.number_of_cells() << " cells expanded c3t3" << std::endl;

    }

    #ifdef CGAL_MESH_3_IO_VERBOSE

      std::cerr << "done.\n";

    #endif

  } // namespace periodic_replicator

} // namespace CGAL

#endif // CGAL_PERIODIC_REPLICATOR_H