Now uses RefFaceList and RefNodeList in MeshNodeBoundary

src/main.cpp

@@ -238,7 +238,6 @@ int main(int argc, char *argv[])
           timer.reset();
           MeshDataType mesh_data(mesh);
 
-
           std::vector<BoundaryConditionHandler> bc_list;
           {
 #warning Should extract boundaries from mesh itself to avoid inconsistencies
@@ -254,7 +253,10 @@ int main(int argc, char *argv[])
                     const RefId& ref = ref_face_list.refId();
                     if (ref == sym_bc_descriptor.boundaryDescriptor()) {
                       std::cout << "Found mesh boundary to impose " << sym_bc_descriptor << '\n';
-                      ref_boundary[ref] = MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list.faceList());
+                      SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+                          = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list));
+                      std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+                      bc_list.push_back(BoundaryConditionHandler(bc));
                     }
                   }
                   break;
@@ -265,48 +267,6 @@ int main(int argc, char *argv[])
                 }
               }
             }
-            for (size_t i_ref_face_list=0; i_ref_face_list<mesh.connectivity().numberOfRefFaceList();
-                 ++i_ref_face_list) {
-              const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
-              const RefId& ref = ref_face_list.refId();
-
-              TinyVector<2> normal(0,0);
-              if ((ref.tagName()== std::string("XMIN")) or (ref.tagName()=="XMAX")) {
-                normal = TinyVector<2>(1,0);
-              } else if ((ref.tagName()=="YMIN") or (ref.tagName()=="YMAX")) {
-                normal = TinyVector<2>(0,1);
-              } else {
-                std::cout << rang::fg::red << "ignoring " << ref << rang::style::reset << '\n';
-                break;
-              }
-
-              std::set<unsigned int> node_id_set;
-              const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes();
-
-              const Kokkos::View<const unsigned int*>& face_id_list = ref_face_list.faceList();
-
-              for (unsigned int l=0; l<face_id_list.extent(0); ++l) {
-                for (unsigned short r=0; r<2; ++r) {
-                  node_id_set.insert(face_nodes(face_id_list[l],r));
-                }
-              }
-
-              Kokkos::View<unsigned int*> node_id_list("node_id_list", node_id_set.size());
-              {
-                int r=0;
-                for (auto node_id : node_id_set) {
-                  node_id_list[r] = node_id;
-                  ++r;
-                }
-              }
-
-              RefNodeList ref_node_list(ref, node_id_list);
-
-              SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
-                  = new SymmetryBoundaryCondition<MeshType::dimension>(ref_node_list, normal);
-              std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
-              bc_list.push_back(BoundaryConditionHandler(bc));
-            }
           }
 
           UnknownsType unknowns(mesh_data);

src/mesh/MeshNodeBoundary.hpp

@@ -4,6 +4,9 @@
 #include <Kokkos_Core.hpp>
 #include <TinyVector.hpp>
 
+#include <RefNodeList.hpp>
+#include <RefFaceList.hpp>
+
 #include <iostream>
 
 template <size_t dimension>
@@ -15,13 +18,19 @@ class MeshNodeBoundary
   MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
   MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
+  const Kokkos::View<const unsigned int*>& nodeList() const
+  {
+    return m_node_list;
+  }
+
   template <typename MeshType>
   MeshNodeBoundary(const MeshType& mesh,
-                   const Kokkos::View<const unsigned int*>& face_list)
+                   const RefFaceList& ref_face_list)
   {
     static_assert(dimension == MeshType::dimension);
     const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells();
 
+    const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
     Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
         if (face_nb_cells[face_list[l]]>1) {
           std::cerr << "internal faces cannot be used to define mesh boundaries\n";
@@ -79,13 +88,18 @@ class MeshFlatNodeBoundary
   template <typename MeshType>
   inline Rd _getOutgoingNormal(const MeshType& mesh);
  public:
+  const Rd& outgoingNormal() const
+  {
+    return m_outgoing_normal;
+  }
+
   MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
   MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
   template <typename MeshType>
   MeshFlatNodeBoundary(const MeshType& mesh,
-                       const Kokkos::View<const unsigned int*>& face_list)
-      : MeshNodeBoundary<dimension>(mesh, face_list),
+                       const RefFaceList& ref_face_list)
+      : MeshNodeBoundary<dimension>(mesh, ref_face_list),
         m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;

src/scheme/BoundaryCondition.hpp

@@ -7,6 +7,7 @@
 #include <Kokkos_Core.hpp>
 
 #include <RefNodeList.hpp>
+#include <MeshNodeBoundary.hpp>
 
 class BoundaryCondition
 {
@@ -84,29 +85,27 @@ public:
   typedef TinyVector<dimension, double> Rd;
 
 private:
-  const RefNodeList m_ref_node_list;
-  // Kokkos::View<unsigned int*> m_node_list;
-  const Rd m_outgoing_normal;
+
+  const MeshFlatNodeBoundary<dimension> m_mesh_flat_node_boundary;
 public:
   const Rd& outgoingNormal() const
   {
-    return m_outgoing_normal;
+    return m_mesh_flat_node_boundary.outgoingNormal();
   }
+
   size_t numberOfNodes() const
   {
-    return m_ref_node_list.nodeList().extent(0);
+    return m_mesh_flat_node_boundary.nodeList().extent(0);
   }
 
   const Kokkos::View<const unsigned int*> nodeList() const
   {
-    return m_ref_node_list.nodeList();
+    return m_mesh_flat_node_boundary.nodeList();
   }
 
-  SymmetryBoundaryCondition(const RefNodeList& ref_node_list,
-                            const Rd& outgoing_normal)
+  SymmetryBoundaryCondition(const MeshFlatNodeBoundary<dimension>& mesh_flat_node_boundary)
     : BoundaryCondition(BoundaryCondition::symmetry),
-      m_ref_node_list(ref_node_list),
-      m_outgoing_normal(outgoing_normal)
+      m_mesh_flat_node_boundary(mesh_flat_node_boundary)
   {
     ;
   }