simple mesh smoother

src/mesh/MeshSmoother.cpp

@@ -167,56 +167,34 @@ class MeshSmootherHandler::MeshSmoother
 
     synchronize(max_delta_xr);
 
-    std::uniform_real_distribution<> distribution(-0.45, 0.45);
-
-    NodeValue<const int> node_numbers = connectivity.nodeNumber();
-    using IdCorrespondance            = std::pair<int, NodeId>;
-    Array<IdCorrespondance> node_numbers_to_node_id{node_numbers.numberOfItems()};
-    parallel_for(
-      node_numbers.numberOfItems(), PUGS_LAMBDA(const NodeId node_id) {
-        node_numbers_to_node_id[node_id] = std::make_pair(node_numbers[node_id], node_id);
-      });
-
-    std::sort(&node_numbers_to_node_id[0], &node_numbers_to_node_id[0] + node_numbers_to_node_id.size(),
-              [](IdCorrespondance a, IdCorrespondance b) { return a.first < b.first; });
-
-    RandomEngine& random_engine = RandomEngine::instance();
-
-    Assert(isSynchronized(random_engine), "seed is not synchronized when entering mesh smoothation");
-
     NodeValue<Rd> shift_r{connectivity};
 
-    int i_node_number = 0;
-    for (size_t i = 0; i < node_numbers_to_node_id.size(); ++i) {
-      const auto [node_number, node_id] = node_numbers_to_node_id[i];
-      while (i_node_number < node_number) {
-        for (size_t j = 0; j < Dimension; ++j) {
-          distribution(random_engine.engine());
-        }
-        ++i_node_number;
-      }
+    parallel_for(
+      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+        const auto& node_cell_list = node_to_cell_matrix[node_id];
 
-      double max_delta = max_delta_xr[node_id];
+        Rd mean_position(zero);
+        size_t number_of_neighbours = 0;
+        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
+          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
 
-      Rd shift;
-      for (size_t i_component = 0; i_component < Dimension; ++i_component) {
-        shift[i_component] = max_delta * distribution(random_engine.engine());
+          const CellId cell_id       = node_cell_list[i_cell];
+          const auto& cell_node_list = cell_to_node_matrix[cell_id];
+          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
+            if (i_node != i_cell_node) {
+              const NodeId cell_node_id = cell_node_list[i_node];
+              mean_position += given_xr[cell_node_id];
+              number_of_neighbours++;
             }
-
-      shift_r[node_id] = shift;
-
-      ++i_node_number;
           }
-
-    const int max_node_number =
-      parallel::allReduceMax(node_numbers_to_node_id[node_numbers_to_node_id.size() - 1].first);
-
-    // Advances random engine to preserve CPU random number generators synchronization
-    for (; i_node_number <= max_node_number; ++i_node_number) {
-      for (size_t j = 0; j < Dimension; ++j) {
-        distribution(random_engine.engine());
         }
+        mean_position    = 1 / number_of_neighbours * mean_position;
+        shift_r[node_id] = mean_position - given_xr[node_id];
+        double nshift    = sqrt(dot(shift_r[node_id], shift_r[node_id]));
+        if (nshift > max_delta_xr[node_id]) {
+          shift_r[node_id] = max_delta_xr[node_id] / nshift * shift_r[node_id];
         }
+      });
 
     this->_applyBC(shift_r);
 
@@ -226,8 +204,6 @@ class MeshSmootherHandler::MeshSmoother
     }
 #endif   // NDEBUG
 
-    Assert(isSynchronized(random_engine), "seed is not synchronized after mesh smoothation");
-
     return shift_r;
   }