diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp index bca86deec4dbb58771cbedc80246e852d4c0fabb..f9a18d44291cd58ab6801c6b515ebab7272f699f 100644 --- a/src/mesh/MeshSmootherEscobar.cpp +++ b/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; @@ -335,8 +346,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar // + trace(transpose(Sigma_gradient[j]) * Sigma_gradient[i]) / Sigma_norm2 // + 0 - 2 * trace(transpose(Sigma) * Sigma_gradient[j]) * trace(transpose(Sigma) * Sigma_gradient[i]) / - (Sigma_norm2 * Sigma_norm2) // - // + (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] // + g[j] * g[i]) * @@ -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;