MeshBoundary now builds node list from face list

MeshFlatboundary uses this list to compute the normal

src/mesh/MeshBoundary.hpp

@@ -10,6 +10,8 @@
 template <size_t dimension>
 class MeshBoundary
 {
+ protected:
+  Kokkos::View<const unsigned int*> m_node_list;
  public:
   MeshBoundary& operator=(const MeshBoundary&) = default;
   MeshBoundary& operator=(MeshBoundary&&) = default;
@@ -27,6 +29,28 @@ class MeshBoundary
           std::exit(1);
         }
       });
+
+    const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes();
+    const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes();
+
+    std::vector<unsigned int> node_ids;
+    // not enough but should reduce significantly the number of resizing
+    node_ids.reserve(dimension*face_list.extent(0));
+    for (size_t l=0; l<face_list.extent(0); ++l) {
+      const size_t face_number = face_list[l];
+      for (size_t r=0; r<face_nb_nodes[face_number]; ++r) {
+        node_ids.push_back(face_nodes(face_number,r));
+      }
+    }
+    std::sort(node_ids.begin(), node_ids.end());
+    auto last = std::unique(node_ids.begin(), node_ids.end());
+    node_ids.resize(std::distance(node_ids.begin(), last));
+
+    Kokkos::View<unsigned int*> node_list("node_list", node_ids.size());
+    Kokkos::parallel_for(node_ids.size(), KOKKOS_LAMBDA(const int& r){
+        node_list[r] = node_ids[r];
+      });
+    m_node_list = node_list;
   }
 
   MeshBoundary() = default;
@@ -45,8 +69,7 @@ class MeshFlatBoundary
   const Rd m_outgoing_normal;
 
   template <typename MeshType>
-  inline Rd _getOutgoingNormal(const MeshType& mesh,
-                               const Kokkos::View<const unsigned int*>& face_list);
+  inline Rd _getOutgoingNormal(const MeshType& mesh);
  public:
   MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default;
   MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default;
@@ -55,7 +78,7 @@ class MeshFlatBoundary
   MeshFlatBoundary(const MeshType& mesh,
                    const Kokkos::View<const unsigned int*>& face_list)
       : MeshBoundary<dimension>(mesh, face_list),
-        m_outgoing_normal(_getOutgoingNormal(mesh, face_list))
+        m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
   }
@@ -71,14 +94,11 @@ template <>
 template <typename MeshType>
 inline TinyVector<2,double>
 MeshFlatBoundary<2>::
-_getOutgoingNormal(const MeshType& mesh,
-                   const Kokkos::View<const unsigned int*>& face_list)
+_getOutgoingNormal(const MeshType& mesh)
 {
   static_assert(MeshType::dimension == 2);
   typedef TinyVector<2,double> R2;
 
-  const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes();
-  const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes();
   const Kokkos::View<const R2*> xr = mesh.xr();
 
   R2 xmin(std::numeric_limits<double>::max(),
@@ -87,15 +107,8 @@ _getOutgoingNormal(const MeshType& mesh,
   R2 xmax(-std::numeric_limits<double>::max(),
           -std::numeric_limits<double>::max());
 
-  if (face_list.extent(0)==0) {
-    std::cerr << "face_list is empty! Cannot compute normal!\n";
-    std::exit(1);
-  }
-
-  for (size_t l=0; l<face_list.extent(0); ++l) {
-    const unsigned int face_number = face_list[l];
-    for (int r=0; r<face_nb_nodes(face_number); ++r) {
-      const R2& x = xr(face_nodes(face_number,r));
+  for (size_t r=0; r<m_node_list.extent(0); ++r) {
+    const R2& x = xr[m_node_list[r]];
     if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
       xmin = x;
     }
@@ -103,7 +116,6 @@ _getOutgoingNormal(const MeshType& mesh,
       xmax = x;
     }
   }
-  }
 
   if (xmin == xmax) {
     std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";