Added calculation of normal to flat boundaries in 3d

src/mesh/MeshNodeBoundary.hpp

@@ -211,6 +211,122 @@ _getNormal(const MeshType& mesh)
   return normal;
 }
 
+template <>
+template <typename MeshType>
+inline TinyVector<3,double>
+MeshFlatNodeBoundary<3>::
+_getNormal(const MeshType& mesh)
+{
+  static_assert(MeshType::dimension == 3);
+  typedef TinyVector<3,double> R3;
+
+
+  R3 xmin(std::numeric_limits<double>::max(),
+          std::numeric_limits<double>::max(),
+          std::numeric_limits<double>::max());
+  R3 ymin = xmin;
+  R3 zmin = xmin;;
+
+  R3 xmax =-xmin;
+  R3 ymax = xmax;
+  R3 zmax = xmax;
+
+  const Kokkos::View<const R3*> xr = mesh.xr();
+
+  for (size_t r=0; r<m_node_list.extent(0); ++r) {
+    const R3& x = xr[m_node_list[r]];
+    if (x[0]<xmin[0]) {
+      xmin = x;
+    }
+    if (x[1]<ymin[1]) {
+      ymin = x;
+    }
+    if (x[2]<zmin[2]) {
+      zmin = x;
+    }
+    if (x[0]>xmax[0]) {
+      xmax = x;
+    }
+    if (x[1]>ymax[1]) {
+      ymax = x;
+    }
+    if (x[2]>zmax[2]) {
+      zmax = x;
+    }
+  }
+
+  const R3 u = xmax-xmin;
+  const R3 v = ymax-ymin;
+  const R3 w = zmax-zmin;
+
+  const R3 uv = crossProduct(u,v);
+  const double uv_l2 = (uv, uv);
+
+  R3 normal = uv;
+  double normal_l2 = uv_l2;
+
+  const R3 uw = crossProduct(u,w);
+  const double uw_l2 = (uw,uw);
+
+  if (uw_l2 > uv_l2) {
+    normal = uw;
+    normal_l2 = uw_l2;
+  }
+
+  const R3 vw = crossProduct(v,w);
+  const double vw_l2 = (vw,vw);
+
+  if (vw_l2 > normal_l2) {
+    normal = vw;
+    normal_l2 = vw_l2;
+  }
+
+
+  if (normal_l2 ==0) {
+    std::cerr << "not able to compute normal!\n";
+    std::exit(1);
+  }
+
+  normal *= 1./sqrt(normal_l2);
+
+#warning Add flatness test
+  // this->_checkBoundaryIsFlat(normal, xmin, xmax, mesh);
+
+  return normal;
+}
+
+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;
+  }
+}
+
 template <>
 template <typename MeshType>
 inline TinyVector<2,double>
@@ -245,16 +361,16 @@ _getOutgoingNormal(const MeshType& mesh)
 
 template <>
 template <typename MeshType>
-inline TinyVector<1,double>
-MeshFlatNodeBoundary<1>::
+inline TinyVector<3,double>
+MeshFlatNodeBoundary<3>::
 _getOutgoingNormal(const MeshType& mesh)
 {
-  static_assert(MeshType::dimension == 1);
-  typedef TinyVector<1,double> R;
+  static_assert(MeshType::dimension == 3);
+  typedef TinyVector<3,double> R3;
 
-  const R normal = this->_getNormal(mesh);
+  const R3 normal = this->_getNormal(mesh);
 
-  const Kokkos::View<const R*>& xr = mesh.xr();
+  const Kokkos::View<const R3*>& 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();