Skip to content
Snippets Groups Projects
Commit 0b1d60dc authored by Stéphane Del Pino's avatar Stéphane Del Pino
Browse files

[ci-skip] Progress in debugging convexifier

parent cd11ade9
No related branches found
No related tags found
No related merge requests found
...@@ -460,7 +460,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -460,7 +460,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
void void
_getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const _getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
{ {
constexpr bool multi_connection = true; constexpr bool multi_connection = false;
const ConnectivityType& connectivity = m_given_mesh.connectivity(); const ConnectivityType& connectivity = m_given_mesh.connectivity();
...@@ -479,7 +479,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -479,7 +479,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
constexpr double eps = 1E-14; constexpr double eps = 1E-14;
quality.fill(2); quality.fill(2);
auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double { size_t maxiter = 0;
auto smooth = [=, &maxiter](const NodeId node_id, TinyVector<Dimension>& x) -> double {
const bool print = (node_id == 4); const bool print = (node_id == 4);
auto cell_list = node_to_cell_matrix[node_id]; auto cell_list = node_to_cell_matrix[node_id];
...@@ -513,7 +515,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -513,7 +515,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
std::cout.precision(15); std::cout.precision(15);
for (size_t i_iter = 0; i_iter < 1; ++i_iter) { for (size_t i_iter = 0; i_iter < 5; ++i_iter) {
double sigma_min = std::numeric_limits<double>::max();
delta = sigma_min;
maxiter = std::max(i_iter + 1, maxiter);
if (print) { if (print) {
std::cout << "------------------------------\n"; std::cout << "------------------------------\n";
std::cout << "iter = " << i_iter << '\n'; std::cout << "iter = " << i_iter << '\n';
...@@ -554,7 +560,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -554,7 +560,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
} }
{ {
double current_sigma_min = min(sigma_list); double current_sigma_min = std::min(min(sigma_list), delta);
if constexpr (multi_connection) { if constexpr (multi_connection) {
size_t nb_medium_range_simplices = 0; size_t nb_medium_range_simplices = 0;
...@@ -858,18 +864,25 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -858,18 +864,25 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
trace(mr_Sigma * mr_S_gradient[1])}; trace(mr_Sigma * mr_S_gradient[1])};
const std::array<TinyMatrix<Dimension>, Dimension> // const std::array<TinyMatrix<Dimension>, Dimension> //
mr_Sigma_gradient{TinyMatrix<Dimension>{0, +1, // mr_Sigma_gradient{TinyMatrix<Dimension>{+0, (1 - mr_cos_alpha) / mr_sin_alpha, //
0, -1}, +0, -1},
TinyMatrix<Dimension>{-mr_inv_sin_alpha, -mr_inv_tan_alpha, // TinyMatrix<Dimension>{(mr_cos_alpha - 1) / mr_sin_alpha, +0, //
+mr_inv_sin_alpha, +mr_inv_tan_alpha}}; +1, +0}};
if (print) { if (print) {
std::cout << "=== MC triangle 2 ===\n"; std::cout << "=== MC triangle 2 ===\n";
std::cout << " A = " << mr_A << '\n'; std::cout << " A = " << mr_A << '\n';
std::cout << " W = " << mr_W << '\n'; std::cout << " W = " << mr_W << '\n';
std::cout << " inv(W) = " << mr_inv_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_S = " << mr_S << '\n';
std::cout << " mr_Sigma = " << mr_Sigma << '\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); escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
...@@ -897,19 +910,24 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -897,19 +910,24 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
const double delta_x_norm = l2Norm(delta_x_candidate); const double delta_x_norm = l2Norm(delta_x_candidate);
double factor = 1;
if (delta_x_norm > 0.3 * max_edge_lenght) {
factor = (0.3 * max_edge_lenght / delta_x_norm);
}
x += factor * delta_x_candidate;
if (print) { if (print) {
std::cout << "==============\n"; std::cout << "==============\n";
std::cout << "i_iter = " << i_iter << '\n';
std::cout << "f_gradient = " << f_gradient << '\n'; std::cout << "f_gradient = " << f_gradient << '\n';
std::cout << "f_hessian = " << f_hessian << '\n'; std::cout << "f_hessian = " << f_hessian << '\n';
std::cout << "inv(f_hes) = " << inverse(f_hessian) << '\n'; std::cout << "inv(f_hes) = " << inverse(f_hessian) << '\n';
std::cout << "delta_x_candidate= " << delta_x_candidate << '\n';
}
double factor = 1;
if (delta_x_norm > 0.2 * max_edge_lenght) {
factor *= (0.2 * max_edge_lenght / delta_x_norm);
}
x += factor * delta_x_candidate;
if (print) {
std::cout << "factor = " << factor << '\n'; std::cout << "factor = " << factor << '\n';
std::cout << "delta_x= " << delta_x_candidate << '\n'; std::cout << "delta_x = " << factor * delta_x_candidate << '\n';
} }
if (l2Norm(f_gradient) < 1E-12) { if (l2Norm(f_gradient) < 1E-12) {
break; break;
...@@ -940,6 +958,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -940,6 +958,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
} }
} }
std::cout << "nb smoothed = " << count << '\n'; std::cout << "nb smoothed = " << count << '\n';
std::cout << "maxiter = " << maxiter << '\n';
} else { } else {
throw NotImplementedError("Dimension != 2"); throw NotImplementedError("Dimension != 2");
} }
...@@ -964,7 +983,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar ...@@ -964,7 +983,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
if (nb_non_convex_cells > 0) { if (nb_non_convex_cells > 0) {
m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node); m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
} }
} while ((nb_non_convex_cells > 0) and (i_convexify < 10)); } while ((nb_non_convex_cells > 0) and (i_convexify < 0));
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr); return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
} }
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment