Added first version of mesh boundaries calculation for BC

src/main.cpp

@@ -23,6 +23,8 @@
 
 #include <BoundaryConditionDescriptor.hpp>
 
+#include <MeshBoundary.hpp>
+
 #include <GmshReader.hpp>
 
 #include <CLI/CLI.hpp>
@@ -239,12 +241,22 @@ int main(int argc, char *argv[])
 
           std::vector<BoundaryConditionHandler> bc_list;
           {
+#warning Should extract boundaries from mesh itself to avoid inconsistencies
+            std::map<RefId, MeshBoundary<MeshType::dimension>> ref_boundary;
             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);
-                  std::cout << sym_bc_descriptor << '\n';
+                  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();
+                    if (ref == sym_bc_descriptor.boundaryDescriptor()) {
+                      std::cout << "Found mesh boundary to impose " << sym_bc_descriptor << '\n';
+                      ref_boundary[ref] = MeshFlatBoundary<MeshType::dimension>(mesh, ref_face_list.faceList());
+                    }
+                  }
                   break;
                 }
                 default: {

src/mesh/Connectivity2D.hpp

@@ -266,6 +266,11 @@ private:
     return m_face_nb_cells;
   }
 
+  const Kokkos::View<const unsigned short*> faceNbNodes() const
+  {
+    return m_face_nb_nodes;
+  }
+
   const Kokkos::View<const unsigned int**> nodeCells() const
   {
     return m_node_cells;

src/mesh/MeshBoundary.hpp

+#ifndef MESH_BOUNDARY_HPP
+#define MESH_BOUNDARY_HPP
+
+#include <Kokkos_Core.hpp>
+#include <TinyVector.hpp>
+
+#include <iostream>
+
+#warning Must change rewrite that to compute associated node lists
+template <size_t dimension>
+class MeshBoundary
+{
+ public:
+  MeshBoundary& operator=(const MeshBoundary&) = default;
+  MeshBoundary& operator=(MeshBoundary&&) = default;
+
+  template <typename MeshType>
+  MeshBoundary(const MeshType& mesh,
+               const Kokkos::View<const unsigned int*>& face_list)
+  {
+    static_assert(dimension == MeshType::dimension);
+    const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells();
+
+    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";
+          std::exit(1);
+        }
+      });
+  }
+
+  MeshBoundary() = default;
+  MeshBoundary(const MeshBoundary&) = default;
+  MeshBoundary(MeshBoundary&&) = default;
+  virtual ~MeshBoundary() = default;
+};
+
+
+template <size_t dimension>
+class MeshFlatBoundary
+    : public MeshBoundary<dimension>
+{
+  typedef TinyVector<dimension, double> Rd;
+ private:
+  const Rd m_outgoing_normal;
+
+  template <typename MeshType>
+  inline Rd _getOutgoingNormal(const MeshType& mesh,
+                               const Kokkos::View<const unsigned int*>& face_list);
+ public:
+  MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default;
+  MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default;
+
+  template <typename MeshType>
+  MeshFlatBoundary(const MeshType& mesh,
+                   const Kokkos::View<const unsigned int*>& face_list)
+      : MeshBoundary<dimension>(mesh, face_list),
+        m_outgoing_normal(_getOutgoingNormal(mesh, face_list))
+  {
+    ;
+  }
+
+  MeshFlatBoundary() = default;
+  MeshFlatBoundary(const MeshFlatBoundary&) = default;
+  MeshFlatBoundary(MeshFlatBoundary&&) = default;
+  virtual ~MeshFlatBoundary() = default;
+
+};
+
+template <>
+template <typename MeshType>
+inline TinyVector<2,double>
+MeshFlatBoundary<2>::
+_getOutgoingNormal(const MeshType& mesh,
+                   const Kokkos::View<const unsigned int*>& face_list)
+{
+  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(),
+          std::numeric_limits<double>::max());
+
+  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));
+      if ((x[0]<xmin[0]) or ((x[0] == xmin[0]) and (x[1] < xmin[1]))) {
+        xmin = x;
+      }
+      if ((x[0]>xmax[0]) or ((x[0] == xmax[0]) and (x[1] > xmax[1]))) {
+        xmax = x;
+      }
+    }
+  }
+
+  if (xmin == xmax) {
+    std::cerr << "xmin==xmax (" << xmin << "==" << xmax << ") unable to compute normal";
+    std::exit(1);
+  }
+
+  R2 dx = xmax-xmin;
+  dx *= 1./std::sqrt((dx,dx));
+
+  R2 normal(-dx[1], dx[0]);
+
+  std::cout << "xmin=" << xmin << " xmax=" << xmax << " normal=" << normal << '\n';
+
+  return normal;
+}
+
+
+#endif // MESH_BOUNDARY_HPP