[ci-skip] Fix grad(Sigma) calculation

src/mesh/MeshSmootherEscobar.cpp

@@ -267,7 +267,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
 
         double final_f = 0;
 
-        for (size_t i_iter = 0; i_iter < 15; ++i_iter) {
+        for (size_t i_iter = 0; i_iter < 100; ++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];
@@ -311,8 +311,19 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
                                                  trace(Sigma * S_gradient[1])};
 
             const std::array<TinyMatrix<Dimension>, Dimension>   //
-              Sigma_gradient{sigma_gradient[0] * inverse(S) - Sigma * S_gradient[0] * Sigma,
-                             sigma_gradient[1] * inverse(S) - Sigma * S_gradient[1] * Sigma};
+              Sigma_gradient_old{sigma_gradient[0] * inverse(S) - inverse(S) * S_gradient[0] * Sigma,
+                                 sigma_gradient[1] * inverse(S) - inverse(S) * S_gradient[1] * Sigma};
+
+            const std::array<TinyMatrix<Dimension>, Dimension>                                           //
+              Sigma_gradient_new{TinyMatrix<Dimension>{0, 1. / std::sin(alpha - 1. / std::tan(alpha)),   //
+                                                       0, -1},
+                                 TinyMatrix<Dimension>{-1. / std::sin(alpha) + 1. / std::tan(alpha), 0,   //
+                                                       1, 0}};
+            const auto Sigma_gradient = Sigma_gradient_new;
+            std::cout << "Sigma_gradient_old[0] = " << Sigma_gradient_old[0] << '\n';
+            std::cout << "Sigma_gradient_new[0] = " << Sigma_gradient_new[0] << '\n';
+            std::cout << "Sigma_gradient_old[1] = " << Sigma_gradient_old[1] << '\n';
+            std::cout << "Sigma_gradient_new[1] = " << Sigma_gradient_new[1] << '\n';
 
             // TinyVector<Dimension> h_gradient = h / (h - sigma_list[i_cell]) * sigma_gradient;
 
@@ -357,7 +368,14 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
           std::cout << "inv(H)         = " << inverse(f_hessian) << '\n';
           std::cout << "inv(H)*grad(f) = " << inverse(f_hessian) * f_gradient << '\n';
 
-          std::cout << "x = " << x << " -> " << x - inverse(f_hessian) * f_gradient << '\n';
+          std::cout << rang::fgB::yellow << "x = " << x << " -> " << x - inverse(f_hessian) * f_gradient
+                    << rang::fg::reset << '\n';
+
+          std::cout << rang::fgB::green << i_iter << ": l2Norm(f_gradient) = " << l2Norm(f_gradient) << rang::fg::reset
+                    << '\n';
+          if (l2Norm(f_gradient) < 1E-6) {
+            break;
+          }
 
           x -= inverse(f_hessian) * f_gradient;