Normal to flat boundaries are now computed

src/main.cpp

@@ -242,7 +242,7 @@ 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;
+            std::map<RefId, MeshNodeBoundary<MeshType::dimension>> ref_boundary;
             for (const auto& bc_descriptor : bc_descriptor_list) {
               switch (bc_descriptor->type()) {
                 case BoundaryConditionDescriptor::Type::symmetry: {
@@ -254,7 +254,7 @@ 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] = MeshFlatBoundary<MeshType::dimension>(mesh, ref_face_list.faceList());
+                      ref_boundary[ref] = MeshFlatNodeBoundary<MeshType::dimension>(mesh, ref_face_list.faceList());
                     }
                   }
                   break;

src/mesh/MeshBoundary.hpp

-#ifndef MESH_BOUNDARY_HPP
-#define MESH_BOUNDARY_HPP
+#ifndef MESH_NODE_BOUNDARY_HPP
+#define MESH_NODE_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
+class MeshNodeBoundary
 {
  protected:
   Kokkos::View<const unsigned int*> m_node_list;
  public:
-  MeshBoundary& operator=(const MeshBoundary&) = default;
-  MeshBoundary& operator=(MeshBoundary&&) = default;
+  MeshNodeBoundary& operator=(const MeshNodeBoundary&) = default;
+  MeshNodeBoundary& operator=(MeshNodeBoundary&&) = default;
 
   template <typename MeshType>
-  MeshBoundary(const MeshType& mesh,
+  MeshNodeBoundary(const MeshType& mesh,
                    const Kokkos::View<const unsigned int*>& face_list)
   {
     static_assert(dimension == MeshType::dimension);
@@ -53,48 +52,80 @@ class MeshBoundary
     m_node_list = node_list;
   }
 
-  MeshBoundary() = default;
-  MeshBoundary(const MeshBoundary&) = default;
-  MeshBoundary(MeshBoundary&&) = default;
-  virtual ~MeshBoundary() = default;
+  MeshNodeBoundary() = default;
+  MeshNodeBoundary(const MeshNodeBoundary&) = default;
+  MeshNodeBoundary(MeshNodeBoundary&&) = default;
+  virtual ~MeshNodeBoundary() = default;
 };
 
 
 template <size_t dimension>
-class MeshFlatBoundary
-    : public MeshBoundary<dimension>
+class MeshFlatNodeBoundary
+    : public MeshNodeBoundary<dimension>
 {
   typedef TinyVector<dimension, double> Rd;
  private:
   const Rd m_outgoing_normal;
 
+  template <typename MeshType>
+  inline Rd _getNormal(const MeshType& mesh);
+
+  template <typename MeshType>
+  inline void _checkBoundaryIsFlat(const TinyVector<2,double>& normal,
+                                   const TinyVector<2,double>& xmin,
+                                   const TinyVector<2,double>& xmax,
+                                   const MeshType& mesh) const;
+
   template <typename MeshType>
   inline Rd _getOutgoingNormal(const MeshType& mesh);
  public:
-  MeshFlatBoundary& operator=(const MeshFlatBoundary&) = default;
-  MeshFlatBoundary& operator=(MeshFlatBoundary&&) = default;
+  MeshFlatNodeBoundary& operator=(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary& operator=(MeshFlatNodeBoundary&&) = default;
 
   template <typename MeshType>
-  MeshFlatBoundary(const MeshType& mesh,
+  MeshFlatNodeBoundary(const MeshType& mesh,
                        const Kokkos::View<const unsigned int*>& face_list)
-      : MeshBoundary<dimension>(mesh, face_list),
+      : MeshNodeBoundary<dimension>(mesh, face_list),
         m_outgoing_normal(_getOutgoingNormal(mesh))
   {
     ;
   }
 
-  MeshFlatBoundary() = default;
-  MeshFlatBoundary(const MeshFlatBoundary&) = default;
-  MeshFlatBoundary(MeshFlatBoundary&&) = default;
-  virtual ~MeshFlatBoundary() = default;
-
+  MeshFlatNodeBoundary() = default;
+  MeshFlatNodeBoundary(const MeshFlatNodeBoundary&) = default;
+  MeshFlatNodeBoundary(MeshFlatNodeBoundary&&) = default;
+  virtual ~MeshFlatNodeBoundary() = default;
 };
 
+template<>
+template <typename MeshType>
+void MeshFlatNodeBoundary<2>::
+_checkBoundaryIsFlat(const TinyVector<2,double>& normal,
+                     const TinyVector<2,double>& xmin,
+                     const TinyVector<2,double>& xmax,
+                     const MeshType& mesh) const
+{
+  static_assert(MeshType::dimension == 2);
+  typedef TinyVector<2,double> R2;
+  const R2 origin = 0.5*(xmin+xmax);
+  const double length = l2Norm(xmax-xmin);
+
+  const Kokkos::View<const R2*> xr = mesh.xr();
+
+  Kokkos::parallel_for(m_node_list.extent(0), KOKKOS_LAMBDA(const size_t& r) {
+      const R2& x = xr[m_node_list[r]];
+      if ((x-origin,normal)>1E-13*length) {
+        std::cerr << "this FlatBoundary is not flat!\n";
+        std::exit(1);
+      }
+    });
+}
+
 template <>
 template <typename MeshType>
 inline TinyVector<2,double>
-MeshFlatBoundary<2>::
-_getOutgoingNormal(const MeshType& mesh)
+MeshFlatNodeBoundary<2>::
+_getNormal(const MeshType& mesh)
 {
   static_assert(MeshType::dimension == 2);
   typedef TinyVector<2,double> R2;
@@ -123,14 +154,46 @@ _getOutgoingNormal(const MeshType& mesh)
   }
 
   R2 dx = xmax-xmin;
-  dx *= 1./std::sqrt((dx,dx));
+  dx *= 1./l2Norm(dx);
 
   R2 normal(-dx[1], dx[0]);
 
-  std::cout << "xmin=" << xmin << " xmax=" << xmax << " normal=" << normal << '\n';
+  this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+
+  return normal;
+}
+
+template <>
+template <typename MeshType>
+inline TinyVector<2,double>
+MeshFlatNodeBoundary<2>::
+_getOutgoingNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::dimension == 2);
+  typedef TinyVector<2,double> R2;
+
+  const R2 normal = this->_getNormal(mesh);
 
+  const Kokkos::View<const R2*>& 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_BOUNDARY_HPP
+#endif // MESH_NODE_BOUNDARY_HPP