Begining of quad mesh support

src/mesh/Connectivity2D.hpp

@@ -131,24 +131,24 @@ public:
     // Computes inefficiently node->cells connectivity [Version 0]
     
     std::multimap<unsigned int, unsigned int> node_cells_map;
+    using namespace Kokkos::Experimental;
+    Kokkos::parallel_reduce(m_number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
+	const size_t n = m_cell_nb_nodes[j];
+	if (n > nb_max) nb_max = n; 
+      }, Max<size_t>(m_max_nb_node_per_cell));
     
-#warning ONLY WORKS FOR TRIANGLES
-    // ------------------------------>>>>
-    m_max_nb_node_per_cell = 3;
     for (unsigned int j=0; j<m_number_of_cells; ++j) {
-      for (unsigned int r=0; r<3; ++r) {
+      for (unsigned int r=0; r<m_cell_nb_nodes[j]; ++r) {
 	node_cells_map.insert(std::make_pair(cell_nodes(j,r),j));
       }
     }
-    // <<<<------------------------------
 
     std::vector<unsigned int> node_ids;
     node_ids.reserve(node_cells_map.size());
     for (const auto& node_cell: node_cells_map) {
       node_ids.push_back(node_cell.first);
     }
-    std::vector<unsigned int>::iterator last_unique
-      = std::unique(node_ids.begin(), node_ids.end());
+    auto last_unique = std::unique(node_ids.begin(), node_ids.end());
     node_ids.resize(std::distance(node_ids.begin(), last_unique));
 
     m_number_of_nodes = node_ids.size();
@@ -156,6 +156,9 @@ public:
     std::cout << "node_ids.size()=" << node_ids.size() << '\n';
     if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) {
       std::cerr << "sparse node numerotation NIY\n";
+      for (int i=0; i<node_ids.size(); ++i) {
+	std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
+      }
       std::exit(0);
     }
 

src/mesh/GmshReader.cpp

@@ -763,24 +763,32 @@ GmshReader::__proceedData()
     if ((dimension1_mask, elementNumber)>0) {
       // throw ErrorHandler(__FILE__,__LINE__,"edges are not treated", ErrorHandler::normal);
     }
-    if (__quadrangles.extent(0)>0) {
-      std::cerr << "quadrangles NYI\n";
-      std::exit(0);
-    }
-
     const size_t nb_cells = (dimension2_mask, elementNumber);
 
     std::cout << "nb_cells=" << nb_cells << '\n';
     const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, 3);
-    for (int j=0; j<nb_cells; ++j) {
+    const size_t nb_triangles = __triangles.extent(0);
+    for (size_t j=0; j<nb_triangles; ++j) {
       cell_nodes(j,0) = __triangles[j][0];
       cell_nodes(j,1) = __triangles[j][1];
       cell_nodes(j,2) = __triangles[j][2];
     }
+    
+    const size_t nb_quadrangles = __quadrangles.extent(0);
+    for (size_t j=0; j<nb_quadrangles; ++j) {
+      const size_t jq = j+nb_triangles;
+      cell_nodes(jq,0) = __quadrangles[j][0];
+      cell_nodes(jq,1) = __quadrangles[j][1];
+      cell_nodes(jq,2) = __quadrangles[j][2];
+      cell_nodes(jq,3) = __quadrangles[j][3];
+    }
     const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
-    for (int j=0; j<nb_cells; ++j) {
+    for (int j=0; j<nb_triangles; ++j) {
       cell_nb_nodes[j] = 3;
     }
+    for (int j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
+      cell_nb_nodes[j] = 4;
+    }
     m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
     Connectivity2D& connectivity = *m_connectivity;
     typedef Mesh<Connectivity2D> MeshType;