From cd11ade9ed089f64caa7cf76b9f7a66f1cf53150 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?St=C3=A9phane=20Del=20Pino?= <stephane.delpino44@gmail.com>
Date: Fri, 28 Jul 2023 19:24:07 +0200
Subject: [PATCH] Begin EMC algorithms [ci-skip]
Debugging still in progress
---
src/mesh/MeshSmootherEscobar.cpp | 395 ++++++++++++++++++++++++++-----
1 file changed, 333 insertions(+), 62 deletions(-)
diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index 412f459b1..eb54ae8fd 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -460,6 +460,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
void
_getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
{
+ constexpr bool multi_connection = true;
+
const ConnectivityType& connectivity = m_given_mesh.connectivity();
auto is_boundary_node = connectivity.isBoundaryNode();
@@ -478,6 +480,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
quality.fill(2);
auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
+ const bool print = (node_id == 4);
+
auto cell_list = node_to_cell_matrix[node_id];
auto node_number_in_cell = node_number_in_their_cells[node_id];
@@ -492,42 +496,30 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const TinyMatrix<Dimension> inv_W = inverse(W);
- std::array<TinyMatrix<Dimension>, Dimension> S_gradient =
- {TinyMatrix<Dimension>{-1, -inv_sin_alpha + inv_tan_alpha, //
- +0, +0}, //
- TinyMatrix<Dimension>{+0, +0, //
- -1, -inv_sin_alpha + inv_tan_alpha}};
+ std::array<TinyMatrix<Dimension>, Dimension> //
+ S_gradient = {TinyMatrix<Dimension>{-1, -1, //
+ +0, +0},
+ TinyMatrix<Dimension>{+inv_tan_alpha, +inv_tan_alpha, //
+ -inv_sin_alpha, -inv_sin_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];
- 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 sigma_min = min(sigma_list);
- double delta =
- (sigma_min < eps) ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min)) : 0;
+ double sigma_min = std::numeric_limits<double>::max();
+ double delta = 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 < 200; ++i_iter) {
+ std::cout.precision(15);
+
+ for (size_t i_iter = 0; i_iter < 1; ++i_iter) {
+ if (print) {
+ std::cout << "------------------------------\n";
+ std::cout << "iter = " << i_iter << '\n';
+ std::cout << "x = " << x << '\n';
+ }
+
double f = 0;
TinyVector<Dimension> f_gradient = zero;
TinyMatrix<Dimension> f_hessian = zero;
@@ -545,6 +537,15 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
a[1] - x[1], b[1] - x[1]};
S_list[i_cell] = A * inv_W;
+
+ if (print) {
+ std::cout << "--- S --- cell = " << i_cell << " (" << node_to_cell_matrix[node_id][i_cell] << ")\n";
+ std::cout << "a = " << a << '\n';
+ std::cout << "b = " << b << '\n';
+ std::cout << "x = " << x << '\n';
+ std::cout << "A = " << A << '\n';
+ std::cout << "S = " << S_list[i_cell] << '\n';
+ }
}
SmallArray<double> sigma_list(S_list.size());
@@ -554,18 +555,102 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
{
double current_sigma_min = min(sigma_list);
+
+ if constexpr (multi_connection) {
+ size_t nb_medium_range_simplices = 0;
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ nb_medium_range_simplices += 2 * (cell_to_node_matrix[cell_id].size() - 3);
+ }
+
+ if (nb_medium_range_simplices > 0) {
+ SmallArray<double> medium_range_sigma_list(nb_medium_range_simplices);
+
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[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 double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
+ const double mr_sin_alpha = std::sin(mr_alpha);
+ const double mr_cos_alpha = std::cos(mr_alpha);
+
+ const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha, //
+ 0, mr_sin_alpha};
+
+ const TinyMatrix mr_inv_W = inverse(mr_W);
+
+ if (print) {
+ std::cout << "MC === i_cell = " << i_cell << '\n';
+ }
+
+ { // first simplex
+ 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 + 2) % cell_nb_nodes]];
+
+ const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
+ a[1] - x[1], b[1] - x[1]};
+
+ const TinyMatrix mr_S = mr_A * mr_inv_W;
+ const double mr_sigma = det(mr_S);
+ if (mr_sigma < current_sigma_min) {
+ current_sigma_min = mr_sigma;
+ }
+
+ if (print) {
+ std::cout << "MC triangle 1 ===\n";
+ std::cout << "MC a = " << a << '\n';
+ std::cout << "MC b = " << b << '\n';
+ std::cout << "MC x = " << x << '\n';
+ }
+ }
+
+ { // second simplex
+ const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
+ const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+
+ const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
+ a[1] - x[1], b[1] - x[1]};
+
+ const TinyMatrix mr_S = mr_A * mr_inv_W;
+ const double mr_sigma = det(mr_S);
+ if (mr_sigma < current_sigma_min) {
+ current_sigma_min = mr_sigma;
+ }
+
+ if (print) {
+ std::cout << "MC triangle 2 ===\n";
+ std::cout << "MC a = " << a << '\n';
+ std::cout << "MC b = " << b << '\n';
+ std::cout << "MC x = " << x << '\n';
+ }
+ }
+ }
+ }
+ }
+
if (current_sigma_min < sigma_min) {
sigma_min = current_sigma_min;
delta = (sigma_min < eps)
? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
: 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];
+ if (print) {
+ std::cout << " ==== \n";
+ std::cout << "delta = " << delta << '\n';
+ }
+ }
+ auto escobarSSigma = [&frobenius, delta, //
+ //
+ &print, &node_to_cell_matrix,
+ &node_id](const double sigma, const TinyMatrix<Dimension>& S,
+ const TinyVector<Dimension>& sigma_gradient,
+ const std::array<TinyMatrix<Dimension>, Dimension>& S_gradient,
+ const std::array<TinyMatrix<Dimension>, Dimension>& Sigma_gradient, double& f,
+ TinyVector<Dimension>& f_gradient, TinyMatrix<Dimension>& f_hessian) {
const TinyMatrix<Dimension> Sigma = [&S] {
TinyMatrix<Dimension> transposed_comS;
for (size_t i = 0; i < Dimension; ++i) {
@@ -585,15 +670,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double f_jr = S_norm * Sigma_norm / h;
- TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
- trace(Sigma * S_gradient[1])};
-
- const std::array<TinyMatrix<Dimension>, Dimension> //
- Sigma_gradient{TinyMatrix<Dimension>{0, inv_sin_alpha - inv_tan_alpha, //
- 0, -1},
- TinyMatrix<Dimension>{-inv_sin_alpha + inv_tan_alpha, 0, //
- 1, 0}};
-
TinyVector<Dimension> g{trace(transpose(S) * S_gradient[0]) / S_norm2 //
+ trace(transpose(Sigma) * Sigma_gradient[0]) / Sigma_norm2 //
- trace(Sigma * S_gradient[0]) / (h - sigma),
@@ -603,7 +679,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
- trace(Sigma * S_gradient[1]) / (h - sigma)};
const TinyVector<Dimension> f_jr_gradient = f_jr * g;
- TinyMatrix<Dimension> f_jr_hessian = zero;
+
+ TinyMatrix<Dimension> f_jr_hessian = zero;
for (size_t i = 0; i < Dimension; ++i) {
for (size_t j = 0; j < Dimension; ++j) {
f_jr_hessian(i, j) = //
@@ -615,20 +692,189 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
- 2 * trace(transpose(Sigma) * Sigma_gradient[j]) * trace(transpose(Sigma) * Sigma_gradient[i]) /
(Sigma_norm2 * Sigma_norm2) //
- - 2 * trace(Sigma_gradient[j] * S_gradient[i]) / (h - sigma) //
- + 2 * trace(Sigma * S_gradient[i]) * sigma / (std::pow(h - sigma, 3)) * sigma_gradient[j] //
+ - trace(Sigma_gradient[j] * S_gradient[i]) / (h - sigma) //
+ + sigma / (std::pow(h - sigma, 3)) * sigma_gradient[i] * sigma_gradient[j] //
+ g[j] * g[i]) *
f_jr;
}
}
+ if (print) {
+ std::cout << "S = " << S << '\n';
+ std::cout << "Sigma = " << Sigma << '\n';
+
+ std::cout << "normS = " << S_norm << '\n';
+ std::cout << "normSigma = " << Sigma_norm << '\n';
+
+ std::cout << "sigma = " << sigma << '\n';
+ std::cout << "f_jr = " << f_jr << '\n';
+ std::cout << "g = " << g << '\n';
+ std::cout << "grad_f_jr = " << f_jr_gradient << '\n';
+ std::cout << "hess_f_jr = " << f_jr_hessian << '\n';
+ }
+
f += f_jr;
f_gradient += f_jr_gradient;
f_hessian += f_jr_hessian;
+ };
+
+ auto escobar = escobarSSigma;
+
+ for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
+ if (print) {
+ std::cout << "------ cell = " << i_cell << '\n';
+ }
+
+ const double sigma = sigma_list[i_cell];
+ const TinyMatrix<Dimension> S = S_list[i_cell];
+
+ const TinyMatrix<Dimension> Sigma = [&S] {
+ TinyMatrix<Dimension> transposed_comS;
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ transposed_comS(i, j) = cofactor(S, j, i);
+ }
+ }
+ return transposed_comS;
+ }();
+
+ TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]), //
+ trace(Sigma * S_gradient[1])};
+
+ const std::array<TinyMatrix<Dimension>, Dimension> //
+ Sigma_gradient{TinyMatrix<Dimension>{+0, +1, //
+ +0, -1},
+ TinyMatrix<Dimension>{-inv_sin_alpha, -inv_tan_alpha, //
+ +inv_sin_alpha, +inv_tan_alpha}};
+
+ escobar(sigma, S, sigma_gradient, S_gradient, Sigma_gradient, f, f_gradient, f_hessian);
}
- if (l2Norm(f_gradient) < 1E-12) {
- break;
+ if constexpr (multi_connection) {
+ for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
+ const CellId cell_id = cell_list[i_cell];
+ auto cell_node_list = cell_to_node_matrix[cell_list[i_cell]];
+ const size_t cell_nb_nodes = cell_node_list.size();
+
+ if (print) {
+ std::cout << "=== MC cell " << i_cell << " ===\n";
+ }
+
+ if (cell_nb_nodes <= 3) {
+ continue;
+ }
+
+ const size_t i_cell_node = node_number_in_cell[i_cell];
+
+ const double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
+ if (print) {
+ std::cout << "=== MC mr_alpha " << mr_alpha << " ===\n";
+ }
+ const double mr_sin_alpha = std::sin(mr_alpha);
+ const double mr_cos_alpha = std::cos(mr_alpha);
+ const double mr_inv_sin_alpha = 1. / mr_sin_alpha;
+ const double mr_inv_tan_alpha = 1. / std::tan(mr_alpha);
+
+ const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha, //
+ 0, mr_sin_alpha};
+
+ const TinyMatrix mr_inv_W = inverse(mr_W);
+
+ std::array<TinyMatrix<Dimension>, Dimension> //
+ mr_S_gradient = {TinyMatrix<Dimension>{-1, (mr_cos_alpha - 1) / mr_sin_alpha, //
+ +0, +0},
+ TinyMatrix<Dimension>{+0, +0, //
+ -1, (mr_cos_alpha - 1) / mr_sin_alpha}};
+
+ { // first simplex
+ 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 + 2) % cell_nb_nodes]];
+
+ const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
+ a[1] - x[1], b[1] - x[1]};
+
+ const TinyMatrix mr_S = mr_A * mr_inv_W;
+ const double mr_sigma = det(mr_S);
+
+ const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
+ TinyMatrix<Dimension> transposed_comS;
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ transposed_comS(i, j) = cofactor(mr_S, j, i);
+ }
+ }
+ return transposed_comS;
+ }();
+
+ TinyVector<Dimension> mr_sigma_gradient{trace(mr_Sigma * mr_S_gradient[0]), //
+ trace(mr_Sigma * mr_S_gradient[1])};
+
+ const std::array<TinyMatrix<Dimension>, Dimension> //
+ mr_Sigma_gradient{TinyMatrix<Dimension>{+0, (1 - mr_cos_alpha) / mr_sin_alpha, //
+ +0, -1},
+ TinyMatrix<Dimension>{(mr_cos_alpha - 1) / mr_sin_alpha, +0, //
+ +1, +0}};
+
+ if (print) {
+ std::cout << "=== MC triangle 1 ===\n";
+ std::cout << " A = " << mr_A << '\n';
+ std::cout << " W = " << mr_W << '\n';
+ std::cout << " inv(W) = " << mr_inv_W << '\n';
+
+ std::cout << " A * inv(W) = " << mr_A * mr_inv_W << '\n';
+
+ std::cout << " mr_S = " << mr_S << '\n';
+ std::cout << " mr_Sigma = " << mr_Sigma << '\n';
+
+ std::cout << " mr_S_gradient = " << mr_S_gradient[0] << " || " << mr_S_gradient[1] << '\n';
+ std::cout << " mr_Sigma_gradient = " << mr_Sigma_gradient[0] << " || " << mr_Sigma_gradient[1]
+ << '\n';
+ }
+
+ escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
+ }
+
+ { // second simplex
+ const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
+ const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+
+ const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0], //
+ a[1] - x[1], b[1] - x[1]};
+
+ const TinyMatrix mr_S = mr_A * mr_inv_W;
+ const double mr_sigma = det(mr_S);
+
+ const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
+ TinyMatrix<Dimension> transposed_comS;
+ for (size_t i = 0; i < Dimension; ++i) {
+ for (size_t j = 0; j < Dimension; ++j) {
+ transposed_comS(i, j) = cofactor(mr_S, j, i);
+ }
+ }
+ return transposed_comS;
+ }();
+
+ TinyVector<Dimension> mr_sigma_gradient{trace(mr_Sigma * mr_S_gradient[0]), //
+ trace(mr_Sigma * mr_S_gradient[1])};
+
+ const std::array<TinyMatrix<Dimension>, Dimension> //
+ mr_Sigma_gradient{TinyMatrix<Dimension>{0, +1, //
+ 0, -1},
+ TinyMatrix<Dimension>{-mr_inv_sin_alpha, -mr_inv_tan_alpha, //
+ +mr_inv_sin_alpha, +mr_inv_tan_alpha}};
+
+ if (print) {
+ std::cout << "=== MC triangle 2 ===\n";
+ std::cout << " A = " << mr_A << '\n';
+ std::cout << " W = " << mr_W << '\n';
+ std::cout << " inv(W) = " << mr_inv_W << '\n';
+ std::cout << " mr_S = " << mr_S << '\n';
+ std::cout << " mr_Sigma = " << mr_Sigma << '\n';
+ }
+
+ escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
+ }
+ }
}
const double max_edge_lenght = [=] {
@@ -651,13 +897,25 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double delta_x_norm = l2Norm(delta_x_candidate);
+ double factor = 1;
if (delta_x_norm > 0.3 * max_edge_lenght) {
- x += (0.3 * max_edge_lenght / delta_x_norm) * delta_x_candidate;
- } else {
- x += delta_x_candidate;
+ factor = (0.3 * max_edge_lenght / delta_x_norm);
+ }
+ x += factor * delta_x_candidate;
+
+ if (print) {
+ std::cout << "==============\n";
+ std::cout << "f_gradient = " << f_gradient << '\n';
+ std::cout << "f_hessian = " << f_hessian << '\n';
+ std::cout << "inv(f_hes) = " << inverse(f_hessian) << '\n';
+ std::cout << "factor = " << factor << '\n';
+ std::cout << "delta_x= " << delta_x_candidate << '\n';
+ }
+ if (l2Norm(f_gradient) < 1E-12) {
+ break;
}
- if (delta_x_norm < 1E-2 * max_edge_lenght) {
+ if (delta_x_norm < 1E-8 * max_edge_lenght) {
break;
}
@@ -667,12 +925,21 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return final_f;
};
- parallel_for(
- connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
- if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
- quality[node_id] = smooth(node_id, new_xr[node_id]);
- }
- });
+ // parallel_for(
+ // connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
+ // if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
+ // quality[node_id] = smooth(node_id, new_xr[node_id]);
+ // }
+ // });
+
+ size_t count = 0;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
+ quality[node_id] = smooth(node_id, new_xr[node_id]);
+ count++;
+ }
+ }
+ std::cout << "nb smoothed = " << count << '\n';
} else {
throw NotImplementedError("Dimension != 2");
}
@@ -690,10 +957,14 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
std::cout << "- convexfication step " << i_convexify++ << '\n';
NodeValue<Rd> new_xr = copy(xr);
this->_getNewPositions(xr, new_xr);
- xr = new_xr;
- nb_non_convex_cells = sum(this->_getNonConvexCellTag(xr));
+ xr = new_xr;
+ auto non_convex_cell_tag = this->_getNonConvexCellTag(xr);
+ nb_non_convex_cells = sum(non_convex_cell_tag);
std::cout << " remaining non convex cells: " << nb_non_convex_cells << '\n';
- } while ((nb_non_convex_cells > 0) and (i_convexify < 120));
+ if (nb_non_convex_cells > 0) {
+ m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
+ }
+ } while ((nb_non_convex_cells > 0) and (i_convexify < 10));
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}
@@ -784,9 +1055,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
- NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
- is_displaced.fill(false);
-
CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
is_zone_cell.fill(false);
@@ -809,6 +1077,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
return m_p_given_mesh;
}
+ NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+ is_displaced.fill(false);
+
for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
auto characteristic_function =
discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
--
GitLab