Plugged symmetry boundary condition descriptor for the 1d case

src/main.cpp

@@ -131,6 +131,17 @@ int main(int argc, char *argv[])
 
       switch (p_mesh->meshDimension()) {
         case 1: {
+          std::vector<std::string> sym_boundary_name_list = {"XMIN", "XMAX"};
+          std::vector<std::shared_ptr<BoundaryConditionDescriptor>> bc_descriptor_list;
+          for (const auto& sym_boundary_name : sym_boundary_name_list){
+            std::shared_ptr<BoundaryDescriptor> boudary_descriptor
+                = std::shared_ptr<BoundaryDescriptor>(new NamedBoundaryDescriptor(sym_boundary_name));
+            SymmetryBoundaryConditionDescriptor* sym_bc_descriptor
+                = new SymmetryBoundaryConditionDescriptor(boudary_descriptor);
+
+            bc_descriptor_list.push_back(std::shared_ptr<BoundaryConditionDescriptor>(sym_bc_descriptor));
+          }
+
           typedef Connectivity1D ConnectivityType;
           typedef Mesh<ConnectivityType> MeshType;
           typedef MeshData<MeshType> MeshDataType;
@@ -143,27 +154,32 @@ int main(int argc, char *argv[])
           MeshDataType mesh_data(mesh);
 
           std::vector<BoundaryConditionHandler> bc_list;
-          // { // quite dirty!
-          // 	for (size_t i_boundary=0; i_boundary<mesh.connectivity().numberOfNodeBoundaries(); ++i_boundary) {
-          // 	  ConnectivityType::NodesBoundary nodes_boundary = mesh.connectivity().nodesBoundary(i_boundary);
-          // 	  unsigned int ref = nodes_boundary.first;
-          // 	  TinyVector<2> normal(0,0);
-          // 	  if ((ref == 5) or (ref == 6)) {
-          // 	    normal = TinyVector<2>(0,1);
-          // 	  } else {
-          // 	    normal = TinyVector<2>(1,0);
-          // 	  }
-          // 	  const Kokkos::View<const unsigned int*> nodes_ids = nodes_boundary.second;
-          // 	  std::vector<unsigned int> node_boundary_vector(nodes_ids.extent(0));
-          // 	  for (size_t r=0; r<nodes_ids.extent(0); ++r) {
-          // 	    node_boundary_vector[r] = nodes_ids[r];
-          // 	  }
-          // 	  SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
-          // 	    = new SymmetryBoundaryCondition<MeshType::dimension>(node_boundary_vector, normal);
-          // 	  std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
-          // 	  bc_list.push_back(BoundaryConditionHandler(bc));
-          // 	}
-          // }
+          {
+            for (const auto& bc_descriptor : bc_descriptor_list) {
+              switch (bc_descriptor->type()) {
+                case BoundaryConditionDescriptor::Type::symmetry: {
+                  const SymmetryBoundaryConditionDescriptor& sym_bc_descriptor
+                      = dynamic_cast<const SymmetryBoundaryConditionDescriptor&>(*bc_descriptor);
+                  for (size_t i_ref_node_list=0; i_ref_node_list<mesh.connectivity().numberOfRefNodeList();
+                       ++i_ref_node_list) {
+                    const RefNodeList& ref_node_list = mesh.connectivity().refNodeList(i_ref_node_list);
+                    const RefId& ref = ref_node_list.refId();
+                    if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                      SymmetryBoundaryCondition<MeshType::dimension>* sym_bc
+                          = new SymmetryBoundaryCondition<MeshType::dimension>(MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_node_list));
+                      std::shared_ptr<SymmetryBoundaryCondition<MeshType::dimension>> bc(sym_bc);
+                      bc_list.push_back(BoundaryConditionHandler(bc));
+                    }
+                  }
+                  break;
+                }
+                default: {
+                  std::cerr << "Unknown BCDescription\n";
+                  std::exit(1);
+                }
+              }
+            }
+          }
 
           UnknownsType unknowns(mesh_data);
 
@@ -212,6 +228,15 @@ int main(int argc, char *argv[])
 
           method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
 
+          { // gnuplot output for density
+            const Kokkos::View<const Rd*> xj   = mesh_data.xj();
+            const Kokkos::View<const double*> rhoj = unknowns.rhoj();
+            std::ofstream fout("rho");
+            for (size_t j=0; j<mesh.numberOfCells(); ++j) {
+              fout << xj[j][0] << ' ' << rhoj[j] << '\n';
+            }
+          }
+
           break;
         }
         case 2: {
@@ -240,8 +265,6 @@ int main(int argc, char *argv[])
 
           std::vector<BoundaryConditionHandler> bc_list;
           {
-#warning Should extract boundaries from mesh itself to avoid inconsistencies
-            std::map<RefId, MeshNodeBoundary<MeshType::dimension>> ref_boundary;
             for (const auto& bc_descriptor : bc_descriptor_list) {
               switch (bc_descriptor->type()) {
                 case BoundaryConditionDescriptor::Type::symmetry: {
@@ -252,7 +275,6 @@ int main(int argc, char *argv[])
                     const RefFaceList& ref_face_list = mesh.connectivity().refFaceList(i_ref_face_list);
                     const RefId& ref = ref_face_list.refId();
                     if (ref == sym_bc_descriptor.boundaryDescriptor()) {
-                      std::cout << "Found mesh boundary to impose " << sym_bc_descriptor << '\n';
                       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);

src/mesh/Connectivity1D.hpp

@@ -7,12 +7,17 @@
 #include <TinyVector.hpp>
 #include <ConnectivityUtils.hpp>
 
+#include <RefId.hpp>
+#include <RefNodeList.hpp>
+
 class Connectivity1D
 {
 public:
   static constexpr size_t dimension = 1;
 
 private:
+  std::vector<RefNodeList> m_ref_node_list;
+
   size_t  m_number_of_nodes;
   const size_t& m_number_of_faces = m_number_of_nodes;
   const size_t  m_number_of_cells;
@@ -61,6 +66,21 @@ private:
   }
 
 public:
+  void addRefNodeList(const RefNodeList& ref_node_list)
+  {
+    m_ref_node_list.push_back(ref_node_list);
+  }
+
+  size_t numberOfRefNodeList() const
+  {
+    return m_ref_node_list.size();
+  }
+
+  const RefNodeList& refNodeList(const size_t& i) const
+  {
+    return m_ref_node_list[i];
+  }
+
   const size_t& numberOfNodes() const
   {
     return m_number_of_nodes;

src/mesh/Connectivity2D.hpp

@@ -19,10 +19,10 @@ class Connectivity2D
 public:
   static constexpr size_t dimension = 2;
 
+private:
   std::vector<RefFaceList> m_ref_face_list;
   std::vector<RefNodeList> m_ref_node_list;
 
-private:
   const size_t  m_number_of_cells;
   size_t m_number_of_faces;
   size_t m_number_of_nodes;

src/mesh/GmshReader.cpp

@@ -841,8 +841,6 @@ GmshReader::__proceedData()
       connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
     }
 
-    const Kokkos::View<const unsigned int**> face_nodes = connectivity.faceNodes();
-
     std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
     for (unsigned int r=0; r<__points.extent(0); ++r) {
       const unsigned int point_number = __points[r];
@@ -885,7 +883,25 @@ GmshReader::__proceedData()
       cell_nb_nodes[j] = 2;
     }
 
-    std::shared_ptr<Connectivity1D> connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes));
+    std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes));
+    Connectivity1D& connectivity = *p_connectivity;
+
+    std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+    for (unsigned int r=0; r<__points.extent(0); ++r) {
+      const unsigned int point_number = __points[r];
+      const unsigned int& ref = __points_ref[r];
+      ref_points_map[ref].push_back(point_number);
+    }
+
+    for (const auto& ref_point_list : ref_points_map) {
+      Kokkos::View<unsigned int*> point_list("point_list", ref_point_list.second.size());
+      for (size_t j=0; j<ref_point_list.second.size(); ++j) {
+        point_list[j]=ref_point_list.second[j];
+      }
+      const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
+      connectivity.addRefNodeList(RefNodeList(physical_ref_id.refId(), point_list));
+    }
+
     typedef Mesh<Connectivity1D> MeshType;
     typedef TinyVector<1, double> Rd;
 
@@ -894,7 +910,7 @@ GmshReader::__proceedData()
       xr[i][0] = __vertices[i][0];
     }
 
-    m_mesh = std::shared_ptr<IMesh>(new MeshType(connectivity, xr));
+    m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
 
   } else {
     std::cerr << "*** using a dimension 0 mesh is forbidden!\n";

src/mesh/MeshNodeBoundary.hpp

@@ -61,6 +61,14 @@ class MeshNodeBoundary
     m_node_list = node_list;
   }
 
+  template <typename MeshType>
+  MeshNodeBoundary(const MeshType& mesh,
+                   const RefNodeList& ref_node_list)
+      : m_node_list(ref_node_list.nodeList())
+  {
+    static_assert(dimension == MeshType::dimension);
+  }
+
   MeshNodeBoundary() = default;
   MeshNodeBoundary(const MeshNodeBoundary&) = default;
   MeshNodeBoundary(MeshNodeBoundary&&) = default;
@@ -105,6 +113,15 @@ class MeshFlatNodeBoundary
     ;
   }
 
+  template <typename MeshType>
+  MeshFlatNodeBoundary(const MeshType& mesh,
+                       const RefNodeList& ref_node_list)
+      : MeshNodeBoundary<dimension>(mesh, ref_node_list),
+        m_outgoing_normal(_getOutgoingNormal(mesh))
+  {
+    ;
+  }
+
   MeshFlatNodeBoundary() = default;
   MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
   MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
@@ -135,6 +152,23 @@ _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
     });
 }
 
+template <>
+template <typename MeshType>
+inline TinyVector<1,double>
+MeshFlatNodeBoundary<1>::
+_getNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::dimension == 1);
+  typedef TinyVector<1,double> R;
+
+  if (m_node_list.extent(0) != 1) {
+    std::cerr << "Node boundaries in 1D require to have exactly 1 node\n";
+    std::exit(1);
+  }
+
+  return R(1);
+}
+
 template <>
 template <typename MeshType>
 inline TinyVector<2,double>
@@ -209,5 +243,37 @@ _getOutgoingNormal(const MeshType& mesh)
   }
 }
 
+template <>
+template <typename MeshType>
+inline TinyVector<1,double>
+MeshFlatNodeBoundary<1>::
+_getOutgoingNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::dimension == 1);
+  typedef TinyVector<1,double> R;
+
+  const R normal = this->_getNormal(mesh);
+
+  const Kokkos::View<const R*>& xr = mesh.xr();
+  const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells();
+  const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes();
+  const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes();
+
+  const size_t r0 = m_node_list[0];
+  const size_t j0 = node_cells(r0,0);
+  double max_height = 0;
+  for (int r=0; r<cell_nb_nodes(j0); ++r) {
+    const double height = (xr(cell_nodes(j0, r))-xr(r0), normal);
+    if (std::abs(height) > std::abs(max_height)) {
+      max_height =  height;
+    }
+  }
+  if (max_height > 0) {
+    return -normal;
+  } else {
+    return normal;
+  }
+}
+
 
 #endif // MESH_NODE_BOUNDARY_HPP