From 02f1e6d711ccad5bb2b4c54102c76fbebd2d38f6 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Mon, 11 Feb 2019 10:31:28 +0100
Subject: [PATCH] Implement *3D* normal calculation in parallel

will require more work to be identical independently of parallelism
---
 src/mesh/MeshNodeBoundary.hpp | 71 +++++++++++++++++++++++++++++------
 1 file changed, 59 insertions(+), 12 deletions(-)

diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 3d470179f..45173f020 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/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,23 +458,37 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R3 normal = this->_getNormal(mesh);
 
-  const NodeValue<const R3>& xr = mesh.xr();
-  const auto& cell_to_node_matrix
-      = mesh.connectivity().cellToNodeMatrix();
+  double max_height = 0;
 
-  const auto& node_to_cell_matrix
-      = mesh.connectivity().nodeToCellMatrix();
+  if (m_node_list.size()>0) {
+    const NodeValue<const R3>& xr = mesh.xr();
+    const auto& cell_to_node_matrix
+        = mesh.connectivity().cellToNodeMatrix();
 
-  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);
+    const auto& node_to_cell_matrix
+        = mesh.connectivity().nodeToCellMatrix();
+
+    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];
+
+    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 {
-- 
GitLab