Fix parallel normal calculation in 1D

src/mesh/MeshNodeBoundary.hpp

@@ -374,6 +374,9 @@ _getOutgoingNormal(const MeshType& mesh)
 
   const R normal = this->_getNormal(mesh);
 
+  double max_height = 0;
+
+  if (m_node_list.size()>0) {
     const NodeValue<const R>& xr = mesh.xr();
     const auto& cell_to_node_matrix
         = mesh.connectivity().cellToNodeMatrix();
@@ -384,13 +387,23 @@ _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 {