diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index 6ef78cc3dd14ab57ee59f7d249979ef123b2abfd..cc36d64587eb52acaa0e7c25cb11f0c3851c036f 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -770,15 +770,91 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
   {
     if constexpr (Dimension == 2) {
       NodeValue<double> node_quality{m_given_mesh.connectivity()};
+
+      auto node_to_cell_matrix        = m_given_mesh.connectivity().nodeToCellMatrix();
+      auto cell_to_node_matrix        = m_given_mesh.connectivity().cellToNodeMatrix();
+      auto node_number_in_their_cells = m_given_mesh.connectivity().nodeLocalNumbersInTheirCells();
+
       for (NodeId node_id = 0; node_id < m_given_mesh.numberOfNodes(); ++node_id) {
         if (not m_is_fixed_node[node_id]) {
-          double sigma_min = std::numeric_limits<double>::max();
-          double delta     = sigma_min;
-          double f;
-          TinyVector<Dimension> f_gradient;
-          TinyMatrix<Dimension> f_hessian;
-          get_f_grad_hess(xr, sigma_min, delta, node_id, xr[node_id], f, f_gradient, f_hessian);
-          node_quality[node_id] = f;
+          auto node_cells  = node_to_cell_matrix[node_id];
+          auto node_number = node_number_in_their_cells[node_id];
+
+          double min_sin     = std::numeric_limits<double>::max();
+          double max_abs_sin = -std::numeric_limits<double>::max();
+
+          for (size_t i_cell = 0; i_cell < node_cells.size(); ++i_cell) {
+            const CellId cell_id = node_cells[i_cell];
+
+            const size_t i_A = node_number[i_cell];
+            const size_t i_B = (i_A + 1) % node_cells.size();
+            const size_t i_C = (i_A + node_cells.size() - 1) % node_cells.size();
+            const size_t i_D = (i_A + 2) % node_cells.size();
+            const size_t i_E = (i_A + node_cells.size() - 2) % node_cells.size();
+
+            const auto cell_nodes = cell_to_node_matrix[cell_id];
+
+            const Rd A = xr[cell_nodes[i_A]];
+            const Rd B = xr[cell_nodes[i_B]];
+            const Rd C = xr[cell_nodes[i_C]];
+            const Rd D = xr[cell_nodes[i_D]];
+            const Rd E = xr[cell_nodes[i_E]];
+
+            const Rd AB = B - A;
+            const Rd AC = C - A;
+            const Rd AD = D - A;
+            const Rd AE = E - A;
+
+            const double l_AB = l2Norm(AB);
+            const double l_AC = l2Norm(AC);
+            const double l_AD = l2Norm(AD);
+            const double l_AE = l2Norm(AE);
+
+            if ((l_AB == 0) or (l_AC == 0) or (l_AD == 0) or (l_AE == 0)) {
+              min_sin     = -1;
+              max_abs_sin = 1;
+              break;
+            }
+
+            const Rd u_AB = 1. / l_AB * AB;
+            const Rd u_AC = 1. / l_AC * AC;
+            const Rd u_AD = 1. / l_AD * AD;
+            const Rd u_AE = 1. / l_AE * AE;
+
+            const double sin_CAB = det(TinyMatrix{u_AB, u_AC});
+            const double sin_DAB = det(TinyMatrix{u_AB, u_AD});
+            const double sin_CAE = det(TinyMatrix{u_AE, u_AC});
+
+            const double abs_sin_CAB = std::abs(sin_CAB);
+            const double abs_sin_DAB = std::abs(sin_DAB);
+            const double abs_sin_CAE = std::abs(sin_CAE);
+
+            if (sin_CAB < min_sin) {
+              min_sin = sin_CAB;
+            }
+
+            if (sin_DAB < min_sin) {
+              min_sin = sin_DAB;
+            }
+
+            if (sin_CAE < min_sin) {
+              min_sin = sin_CAE;
+            }
+
+            if (abs_sin_CAB > max_abs_sin) {
+              max_abs_sin = abs_sin_CAB;
+            }
+
+            if (abs_sin_DAB > max_abs_sin) {
+              max_abs_sin = abs_sin_DAB;
+            }
+
+            if (abs_sin_CAE > max_abs_sin) {
+              max_abs_sin = abs_sin_CAE;
+            }
+          }
+
+          node_quality[node_id] = min_sin / max_abs_sin;
         } else {
           node_quality[node_id] = 1;
         }
@@ -1080,10 +1156,27 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
 
       auto node_quality = this->_getNodeQuality(xr);
 
-      bool need_smooth = max(node_quality) > 50;
+      CellValue<const int> non_convex_cell_tag = [=, this] {
+        CellValue<int> cell_tag(m_given_mesh.connectivity());
+        cell_tag.fill(0);
+
+        auto cell_to_node_matrix = m_given_mesh.connectivity().cellToNodeMatrix();
+
+        for (CellId cell_id = 0; cell_id < m_given_mesh.numberOfCells(); ++cell_id) {
+          for (size_t i_node = 0; i_node < cell_to_node_matrix[cell_id].size(); ++i_node) {
+            const NodeId node_id = cell_to_node_matrix[cell_id][i_node];
+            if (node_quality[node_id] <= 0.15) {
+              cell_tag[cell_id] = 1;
+              break;
+            }
+          }
+        }
+
+        return cell_tag;
+      }();
 
-      auto non_convex_cell_tag = this->_getNonConvexCellTag(xr);
-      nb_non_convex_cells      = sum(non_convex_cell_tag);
+      //   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';
       if (nb_non_convex_cells > 0) {
@@ -1091,7 +1184,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
         this->_getNewPositions(xr, new_xr);
         xr = new_xr;
       }
-    } while ((nb_non_convex_cells > 0) and (i_convexify < 10000));
+    } while ((nb_non_convex_cells > 0) and (i_convexify < 20));
 
     return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
   }