Change maximum displacement calculation

It is now computed according to nodes of the neighboring cells and not
only using the neighboring edges

src/mesh/MeshRandomizer.hpp

@@ -143,29 +143,34 @@ class MeshRandomizer
     const ConnectivityType& connectivity = m_given_mesh.connectivity();
     NodeValue<const Rd> given_xr         = m_given_mesh.xr();
 
-    auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
-    EdgeValue<double> edge_length{connectivity};
-    parallel_for(
-      connectivity.numberOfEdges(), PUGS_LAMBDA(const EdgeId edge_id) {
-        const auto& edge_nodes = edge_to_node_matrix[edge_id];
-        const NodeId node_0    = edge_nodes[0];
-        const NodeId node_1    = edge_nodes[1];
-        const Rd delta_x       = given_xr[node_0] - given_xr[node_1];
-
-        edge_length[edge_id] = std::sqrt(dot(delta_x, delta_x));
-      });
+    auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
+    auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
+    auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
 
-    auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
     NodeValue<double> max_delta_xr{connectivity};
     parallel_for(
       connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        const auto& node_edges = node_to_edge_matrix[node_id];
-        double max_delta       = edge_length[node_edges[0]];
+        const Rd& x0 = given_xr[node_id];
+
+        const auto& node_cell_list = node_to_cell_matrix[node_id];
+        double min_distance_2      = std::numeric_limits<double>::max();
+
+        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);
 
-        for (size_t i_edge = 1; i_edge < node_edges.size(); ++i_edge) {
-          const EdgeId edge_id = node_edges[i_edge];
-          max_delta            = std::min(max_delta, edge_length[edge_id]);
+          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];
+              const Rd delta            = x0 - given_xr[cell_node_id];
+              min_distance_2            = std::min(min_distance_2, dot(delta, delta));
+            }
           }
+        }
+        double max_delta = std::sqrt(min_distance_2);
+
         max_delta_xr[node_id] = max_delta;
       });