diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index eb54ae8fd5e5ba5432e6b0e5bfa88e799d1e4bac..5f608eb2e668d4a7c0d6914d40ca813c589fde69 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -263,8 +263,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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;
+ ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
+ : 0;
}
}
@@ -460,7 +460,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
void
_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();
@@ -479,7 +479,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
constexpr double eps = 1E-14;
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);
auto cell_list = node_to_cell_matrix[node_id];
@@ -513,7 +515,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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) {
std::cout << "------------------------------\n";
std::cout << "iter = " << i_iter << '\n';
@@ -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) {
size_t nb_medium_range_simplices = 0;
@@ -633,8 +639,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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;
+ ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
+ : 0;
}
if (print) {
@@ -857,19 +863,26 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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}};
+ 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 2 ===\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);
@@ -897,19 +910,24 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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) {
std::cout << "==============\n";
+ std::cout << "i_iter = " << i_iter << '\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 << "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 << "delta_x= " << delta_x_candidate << '\n';
+ std::cout << "delta_x = " << factor * delta_x_candidate << '\n';
}
if (l2Norm(f_gradient) < 1E-12) {
break;
@@ -940,6 +958,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
}
}
std::cout << "nb smoothed = " << count << '\n';
+ std::cout << "maxiter = " << maxiter << '\n';
} else {
throw NotImplementedError("Dimension != 2");
}
@@ -964,7 +983,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
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));
+ } while ((nb_non_convex_cells > 0) and (i_convexify < 0));
return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
}