[ci-skip] Add symmetry treatment in 2D (still buggy)

src/mesh/MeshSmootherEscobar.cpp

@@ -99,7 +99,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
         boundary_condition);
     }
 
-#warning treat the axis line case in dimension 3 before this
+    if constexpr (Dimension == 3) {
+      throw NotImplementedError("treat sliding axis node kind");
+    }
 
     synchronize(is_fixed);
 
@@ -164,18 +166,164 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
           if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
             const Rd& n = bc.outgoingNormal();
 
+            const auto node_list = bc.nodeList();
+
             const Rdxd I   = identity;
             const Rdxd nxn = tensorProduct(n, n);
             const Rdxd P   = I - nxn;
 
-            const Array<const NodeId>& node_list = bc.nodeList();
+            const ConnectivityType& connectivity = m_given_mesh.connectivity();
+
+            auto is_boundary_node = connectivity.isBoundaryNode();
+            auto node_is_owned    = connectivity.nodeIsOwned();
+
+            if constexpr (Dimension == 2) {
+              const TinyVector<Dimension> t{-n[1], n[0]};
+
+              auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
+              auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
+              auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
+
+              NodeValue<double> quality{connectivity};
+
+              constexpr double eps = 1E-15;
+              quality.fill(2);
+
+              auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
+                auto cell_list           = node_to_cell_matrix[node_id];
+                auto node_number_in_cell = node_number_in_their_cells[node_id];
+
+                const double alpha = std::acos(-1) / cell_list.size();
+                const TinyMatrix<Dimension> W{1, std::cos(alpha),   //
+                                              0, std::sin(alpha)};
+
+                const TinyMatrix<Dimension> inv_W = inverse(W);
+
+                TinyMatrix<Dimension, Dimension> dtheta_S =
+                  TinyMatrix<Dimension>{t[0], t[0] * (1. / std::sin(alpha) - 1. / std::tan(alpha)),   //
+                                        t[1], t[1] * (1. / std::sin(alpha) - 1. / std::tan(alpha))};
+
+                SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
+                for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+                  const size_t i_cell_node   = node_number_in_cell[i_cell];
+                  const auto cell_node_list  = cell_to_node_matrix[cell_list[i_cell]];
+                  const size_t cell_nb_nodes = cell_node_list.size();
+
+                  const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
+                  const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+
+                  const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0],   //
+                                                a[1] - x[1], b[1] - x[1]};
+
+                  S_list[i_cell] = A * inv_W;
+                }
+
+                SmallArray<double> sigma_list(S_list.size());
+                for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
+                  sigma_list[i_cell] = det(S_list[i_cell]);
+                }
+
+                const double sigma_min = min(sigma_list);
+
+                const double delta =
+                  (sigma_min < eps) ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
+                                    : 0;
+
+                auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
+
+                double final_f = 0;
+
+                for (size_t i_iter = 0; i_iter < 3; ++i_iter) {
+                  SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
+                  for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+                    const size_t i_cell_node   = node_number_in_cell[i_cell];
+                    auto cell_node_list        = cell_to_node_matrix[cell_list[i_cell]];
+                    const size_t cell_nb_nodes = cell_node_list.size();
+
+                    const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
+                    const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+
+                    const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0],   //
+                                                  a[1] - x[1], b[1] - x[1]};
+
+                    S_list[i_cell] = A * inv_W;
+                  }
+
+                  SmallArray<double> sigma_list(S_list.size());
+                  for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
+                    sigma_list[i_cell] = det(S_list[i_cell]);
+                  }
+
+                  double f         = 0;
+                  double dtheta_f  = 0;
+                  double d2theta_f = 0;
+
+                  for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
+                    const double sigma            = sigma_list[i_cell];
+                    const TinyMatrix<Dimension> S = S_list[i_cell];
+
+                    const TinyMatrix<Dimension> Sigma = sigma * inverse(S);
+
+                    const double S_norm      = frobenius(S);
+                    const double Sigma_norm  = frobenius(Sigma);
+                    const double S_norm2     = S_norm * S_norm;
+                    const double Sigma_norm2 = Sigma_norm * Sigma_norm;
+
+                    const double h = sigma + std::sqrt(sigma * sigma + 4 * delta * delta);
+
+                    const double f_jr = S_norm * Sigma_norm / h;
+
+                    const double dtheta_sigma = trace(Sigma * dtheta_S);
+
+                    TinyMatrix<Dimension, Dimension> dtheta_Sigma =
+                      TinyMatrix<Dimension>{t[1] * (1. / std::sin(alpha) - 1. / std::tan(alpha)),
+                                            -t[0] * (1. / std::sin(alpha) - 1. / std::tan(alpha)),   //
+                                            -t[1], t[0]};
+
+                    const double g = trace(transpose(S) * dtheta_S) / S_norm2                 //
+                                     + trace(transpose(Sigma) * dtheta_Sigma) / Sigma_norm2   //
+                                     - dtheta_sigma / (h - sigma);
+
+                    const double dtheta_f_jr = f_jr * g;
+                    const double dtheta2_f_jr =
+                      (trace(transpose(dtheta_S) * dtheta_S) / S_norm2                                               //
+                       - 2 * trace(transpose(S) * dtheta_S) * trace(transpose(S) * dtheta_S) / (S_norm2 * S_norm2)   //
+                                                                                                                     //
+                       + trace(transpose(dtheta_Sigma) * dtheta_Sigma) / Sigma_norm2   // + 0
+                       - 2 * trace(transpose(Sigma) * dtheta_Sigma) * trace(transpose(Sigma) * dtheta_Sigma) /
+                           (Sigma_norm2 * Sigma_norm2)   //
+                       //
+                       - 2 * trace(dtheta_Sigma * dtheta_S) / (h - sigma)                     //
+                       + 2 * sigma / (std::pow(h - sigma, 3)) * dtheta_sigma * dtheta_sigma   //
+                       + g * g) *
+                      f_jr;
+
+                    f += f_jr;
+                    dtheta_f += dtheta_f_jr;
+                    d2theta_f += dtheta2_f_jr;
+                  }
+                  if (std::abs(dtheta_f) < 1E-6) {
+                    break;
+                  }
+                  x += dtheta_f / d2theta_f * t;
+
+                  final_f = f;
+                }
+                return final_f;
+              };
+
               parallel_for(
                 node_list.size(), PUGS_LAMBDA(const size_t i_node) {
                   const NodeId node_id = node_list[i_node];
-
-                new_xr[node_id] = P * new_xr[node_id];
+                  if (not m_is_fixed_node[node_id] and node_is_owned[node_id]) {
+                    smooth(node_id, new_xr[node_id]);
+                  }
                 });
 
+            } else {
+              throw NotImplementedError("Dimension != 2");
+            }
+
           } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
             if constexpr (Dimension > 1) {
               throw NotImplementedError("Escobar: axis boundary conditions");