diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index 736a0e03725e1910605779898e3107ff82575144..c1ff9e69feb146b15716498ec72c28c3c40e4752 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/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,17 +166,163 @@ 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();
- parallel_for(
- node_list.size(), PUGS_LAMBDA(const size_t i_node) {
- const NodeId node_id = node_list[i_node];
+ const ConnectivityType& connectivity = m_given_mesh.connectivity();
- new_xr[node_id] = P * new_xr[node_id];
- });
+ 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];
+ 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) {