Implement *3D* normal calculation in parallel

will require more work to be identical independently of parallelism

src/mesh/MeshNodeBoundary.hpp

@@ -183,7 +183,7 @@ _getNormal(const MeshType&)
     std::exit(1);
   }
 
-  return R(1);
+  return R{1};
 }
 
 template <>
@@ -288,6 +288,39 @@ _getNormal(const MeshType& mesh)
       zmax = x;
     }
   }
+#warning re work this part to avoir parallelism dependance
+  Array<R3> xmin_array = parallel::allGather(xmin);
+  Array<R3> xmax_array = parallel::allGather(xmax);
+  Array<R3> ymin_array = parallel::allGather(ymin);
+  Array<R3> ymax_array = parallel::allGather(ymax);
+  Array<R3> zmin_array = parallel::allGather(zmin);
+  Array<R3> zmax_array = parallel::allGather(zmax);
+
+  for (size_t i=0; i<xmin_array.size(); ++i) {
+    const R3& x = xmin_array[i];
+    if (x[0] < xmin[0]) { xmin = x; }
+  }
+  for (size_t i=0; i<ymin_array.size(); ++i) {
+    const R3& x = ymin_array[i];
+    if (x[1] < ymin[1]) { ymin = x; }
+  }
+  for (size_t i=0; i<zmin_array.size(); ++i) {
+    const R3& x = zmin_array[i];
+    if (x[2] < zmin[2]) { zmin = x; }
+  }
+  for (size_t i=0; i<xmax_array.size(); ++i) {
+    const R3& x = xmax_array[i];
+    if (x[0] > xmax[0]) { xmax = x; }
+  }
+  for (size_t i=0; i<ymax_array.size(); ++i) {
+    const R3& x = ymax_array[i];
+    if (x[1] > ymax[1]) { ymax = x; }
+  }
+  for (size_t i=0; i<zmax_array.size(); ++i) {
+    const R3& x = zmax_array[i];
+    if (x[2] > zmax[2]) { zmax = x; }
+  }
+
 
   const R3 u = xmax-xmin;
   const R3 v = ymax-ymin;
@@ -425,6 +458,9 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R3 normal = this->_getNormal(mesh);
 
+  double max_height = 0;
+
+  if (m_node_list.size()>0) {
     const NodeValue<const R3>& xr = mesh.xr();
     const auto& cell_to_node_matrix
         = mesh.connectivity().cellToNodeMatrix();
@@ -435,13 +471,24 @@ _getOutgoingNormal(const MeshType& mesh)
     const NodeId r0 = m_node_list[0];
     const CellId j0 = node_to_cell_matrix[r0][0];
     const auto& j0_nodes = cell_to_node_matrix[j0];
-  double max_height = 0;
+
     for (size_t r=0; r<j0_nodes.size(); ++r) {
       const double height = (xr[j0_nodes[r]]-xr[r0], normal);
       if (std::abs(height) > std::abs(max_height)) {
         max_height = height;
       }
     }
+
+  }
+
+  Array<double> max_height_array = parallel::allGather(max_height);
+  for (size_t i=0; i<max_height_array.size(); ++i) {
+    const double height = max_height_array[i];
+    if (std::abs(height) > std::abs(max_height)) {
+      max_height = height;
+    }
+  }
+
   if (max_height > 0) {
     return -normal;
   } else {