diff --git a/src/mesh/MeshSmootherEscobar.cpp b/src/mesh/MeshSmootherEscobar.cpp
index f9a18d44291cd58ab6801c6b515ebab7272f699f..736a0e03725e1910605779898e3107ff82575144 100644
--- a/src/mesh/MeshSmootherEscobar.cpp
+++ b/src/mesh/MeshSmootherEscobar.cpp
@@ -40,9 +40,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
   using BoundaryConditionList = std::vector<BoundaryCondition>;
   BoundaryConditionList m_boundary_condition_list;
 
+  NodeValue<const bool> m_is_fixed_node;
+
   BoundaryConditionList
   _getBCList(const MeshType& mesh,
-             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
+             const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
   {
     BoundaryConditionList bc_list;
 
@@ -77,8 +79,83 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
     return bc_list;
   }
 
+  NodeValue<const bool>
+  _getFixedNodes() const
+  {
+    NodeValue<bool> is_fixed{m_given_mesh.connectivity()};
+    is_fixed.fill(false);
+
+    for (auto&& boundary_condition : m_boundary_condition_list) {
+      std::visit(
+        [&](auto&& bc) {
+          using BCType = std::decay_t<decltype(bc)>;
+          if constexpr ((std::is_same_v<BCType, FixedBoundaryCondition>) or
+                        ((Dimension == 2) and std::is_same_v<BCType, AxisBoundaryCondition>)) {
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              node_list.size(), PUGS_LAMBDA(size_t i_node) { is_fixed[node_list[i_node]] = true; });
+          }
+        },
+        boundary_condition);
+    }
+
+#warning treat the axis line case in dimension 3 before this
+
+    synchronize(is_fixed);
+
+    NodeValue<int> bc_number{m_given_mesh.connectivity()};
+    bc_number.fill(-1);
+
+    {
+      int bc_id = 0;
+
+      for (auto&& boundary_condition : m_boundary_condition_list) {
+        std::visit(
+          [&](auto&& bc) {
+            using BCType = std::decay_t<decltype(bc)>;
+            if constexpr (std::is_same_v<BCType, SymmetryBoundaryCondition>) {
+              const Array<const NodeId>& node_list = bc.nodeList();
+              parallel_for(
+                node_list.size(), PUGS_LAMBDA(size_t i_node) {
+                  const NodeId node_id = node_list[i_node];
+                  if (not is_fixed[node_id]) {
+                    if (bc_number[node_id] < 0) {
+                      bc_number[node_id] = bc_id;
+                    } else {
+                      bc_number[node_id] = -2;   // error tag
+                    }
+                  }
+                });
+            }
+          },
+          boundary_condition);
+        bc_id++;
+      }
+
+      if (min(bc_number) < -1) {
+        throw NormalError("Smoothing boundary conditions overlap. One must precise fixed points and smoothing axis");
+      }
+    }
+
+    bool missing_bc  = false;
+    auto is_boundary = m_given_mesh.connectivity().isBoundaryNode();
+    parallel_for(m_given_mesh.numberOfNodes(), [&](NodeId node_id) {
+      if (is_boundary[node_id]) {
+        if (not is_fixed[node_id] and (bc_number[node_id] == -1)) {
+          missing_bc = true;
+        }
+      }
+    });
+
+    if (parallel::allReduceOr(missing_bc)) {
+      throw NormalError("Some boundary conditions are missing (boundary nodes without BC)");
+    }
+
+    return is_fixed;
+  }
+
   void
-  _applyBC(NodeValue<Rd>& shift) const
+  _applyBC(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
   {
     for (auto&& boundary_condition : m_boundary_condition_list) {
       std::visit(
@@ -96,22 +173,12 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
               node_list.size(), PUGS_LAMBDA(const size_t i_node) {
                 const NodeId node_id = node_list[i_node];
 
-                shift[node_id] = P * shift[node_id];
+                new_xr[node_id] = P * new_xr[node_id];
               });
 
           } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
             if constexpr (Dimension > 1) {
-              const Rd& t = bc.direction();
-
-              const Rdxd txt = tensorProduct(t, t);
-
-              const Array<const NodeId>& node_list = bc.nodeList();
-              parallel_for(
-                node_list.size(), PUGS_LAMBDA(const size_t i_node) {
-                  const NodeId node_id = node_list[i_node];
-
-                  shift[node_id] = txt * shift[node_id];
-                });
+              throw NotImplementedError("Escobar: axis boundary conditions");
             } else {
               throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
             }
@@ -121,7 +188,7 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
             parallel_for(
               node_list.size(), PUGS_LAMBDA(const size_t i_node) {
                 const NodeId node_id = node_list[i_node];
-                shift[node_id]       = zero;
+                new_xr[node_id]      = old_xr[node_id];
               });
 
           } else {
@@ -132,94 +199,32 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
     }
   }
 
-  NodeValue<Rd>
-  _getDisplacement() const
-  {
-    const ConnectivityType& connectivity = m_given_mesh.connectivity();
-    NodeValue<const Rd> given_xr         = m_given_mesh.xr();
-
-    auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
-    auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
-    auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
-
-    NodeValue<double> max_delta_xr{connectivity};
-    parallel_for(
-      connectivity.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        const Rd& x0 = given_xr[node_id];
-
-        const auto& node_cell_list = node_to_cell_matrix[node_id];
-        double min_distance_2      = std::numeric_limits<double>::max();
-
-        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
-          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
-
-          const CellId cell_id       = node_cell_list[i_cell];
-          const auto& cell_node_list = cell_to_node_matrix[cell_id];
-
-          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
-            if (i_node != i_cell_node) {
-              const NodeId cell_node_id = cell_node_list[i_node];
-              const Rd delta            = x0 - given_xr[cell_node_id];
-              min_distance_2            = std::min(min_distance_2, dot(delta, delta));
-            }
-          }
-        }
-        double max_delta = std::sqrt(min_distance_2);
-
-        max_delta_xr[node_id] = max_delta;
-      });
-
-    NodeValue<Rd> shift_r{connectivity};
-
-    parallel_for(
-      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-        const auto& node_cell_list = node_to_cell_matrix[node_id];
-        Rd mean_position(zero);
-        size_t number_of_neighbours = 0;
-
-        for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
-          const size_t i_cell_node = node_number_in_their_cells(node_id, i_cell);
-
-          const CellId cell_id       = node_cell_list[i_cell];
-          const auto& cell_node_list = cell_to_node_matrix[cell_id];
-          for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
-            if (i_node != i_cell_node) {
-              const NodeId cell_node_id = cell_node_list[i_node];
-              mean_position += given_xr[cell_node_id];
-              number_of_neighbours++;
-            }
-          }
-        }
-        mean_position    = 1. / number_of_neighbours * mean_position;
-        shift_r[node_id] = mean_position - given_xr[node_id];
-      });
-
-    this->_applyBC(shift_r);
-
-    synchronize(shift_r);
-
-    return shift_r;
-  }
-
  public:
   std::shared_ptr<const ItemValueVariant>
   getQuality() const
   {
-    if constexpr (Dimension == 2) {
-      const ConnectivityType& connectivity = m_given_mesh.connectivity();
-      NodeValue<const Rd> xr               = m_given_mesh.xr();
+    return std::make_shared<ItemValueVariant>(m_given_mesh.xr());
+  }
 
+  void
+  _getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
+  {
+    const ConnectivityType& connectivity = m_given_mesh.connectivity();
+
+    auto is_boundary_node = connectivity.isBoundaryNode();
+    auto node_is_owned    = connectivity.nodeIsOwned();
+
+    if constexpr (Dimension == 2) {
       auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
       auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
       auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
 
-      auto is_boundary_node = connectivity.isBoundaryNode();
       NodeValue<double> quality{connectivity};
 
       constexpr double eps = 1E-15;
       quality.fill(2);
 
-      auto f_inner = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
+      auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
         auto cell_list           = node_to_cell_matrix[node_id];
         auto node_number_in_cell = node_number_in_their_cells[node_id];
 
@@ -241,8 +246,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
           auto cell_node_list        = cell_to_node_matrix[cell_list[i_cell]];
           const size_t cell_nb_nodes = cell_node_list.size();
 
-          const TinyVector a = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
-          const TinyVector b = xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+          const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
+          const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
 
           const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0],   //
                                         a[1] - x[1], b[1] - x[1]};
@@ -262,9 +267,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
 
         auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
 
-        // TinyVector<Dimension> f_gradient = zero;
-        // TinyMatrix<Dimension> f_hessian  = zero;
-
         double final_f = 0;
 
         for (size_t i_iter = 0; i_iter < 100; ++i_iter) {
@@ -274,8 +276,8 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
             auto cell_node_list        = cell_to_node_matrix[cell_list[i_cell]];
             const size_t cell_nb_nodes = cell_node_list.size();
 
-            const TinyVector a = xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
-            const TinyVector b = xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
+            const TinyVector a = old_xr[cell_node_list[(i_cell_node + 1) % cell_nb_nodes]];
+            const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
 
             const TinyMatrix<Dimension> A{a[0] - x[0], b[0] - x[0],   //
                                           a[1] - x[1], b[1] - x[1]};
@@ -310,22 +312,11 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
             TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]),   //
                                                  trace(Sigma * S_gradient[1])};
 
-            const std::array<TinyMatrix<Dimension>, Dimension>   //
-              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;
+            const std::array<TinyMatrix<Dimension>, Dimension>                                       //
+              Sigma_gradient{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}};
 
             TinyVector<Dimension> g{trace(transpose(S) * S_gradient[0]) / S_norm2                   //
                                       + trace(transpose(Sigma) * Sigma_gradient[0]) / Sigma_norm2   //
@@ -360,19 +351,6 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
             f_hessian += f_jr_hessian;
           }
 
-          std::cout << "f = " << f << '\n';
-          std::cout << "grad(f) = " << f_gradient << '\n';
-          std::cout << "hess(f) = " << f_hessian << " | hess(f)^T -hess(f) = " << transpose(f_hessian) - f_hessian
-                    << '\n';
-
-          std::cout << "inv(H)         = " << inverse(f_hessian) << '\n';
-          std::cout << "inv(H)*grad(f) = " << 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;
           }
@@ -386,132 +364,128 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
 
       parallel_for(
         connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
-          // auto cell_list           = node_to_cell_matrix[node_id];
-          // auto node_number_in_cell = node_number_in_their_cells[node_id];
-
-          if (is_boundary_node[node_id]) {
-            quality[node_id] = 1;
-          } else {
-            TinyVector x     = xr[node_id];
-            quality[node_id] = f_inner(node_id, x);
-
-            std::exit(0);
-
-            // TinyMatrix<Dimension> B = identity;
+          if (node_is_owned[node_id] and not is_boundary_node[node_id]) {
+            quality[node_id] = smooth(node_id, new_xr[node_id]);
           }
         });
-
-      return std::make_shared<ItemValueVariant>(quality);
     } else {
       throw NotImplementedError("Dimension != 2");
     }
+
+    this->_applyBC(old_xr, new_xr);
   }
 
   std::shared_ptr<const IMesh>
   getSmoothedMesh() const
   {
     NodeValue<const Rd> given_xr = m_given_mesh.xr();
+    NodeValue<Rd> new_xr         = copy(given_xr);
 
-    NodeValue<Rd> xr = this->_getDisplacement();
-
-    parallel_for(
-      m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) { xr[node_id] += given_xr[node_id]; });
+    this->_getNewPositions(given_xr, new_xr);
 
-    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), new_xr);
   }
 
   std::shared_ptr<const IMesh>
   getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
   {
-    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+    // NodeValue<const Rd> given_xr = m_given_mesh.xr();
 
-    NodeValue<const bool> is_displaced =
-      InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
+    // NodeValue<const bool> is_displaced =
+    //   InterpolateItemValue<bool(const Rd)>::interpolate(function_symbol_id, given_xr);
 
-    NodeValue<Rd> xr = this->_getDisplacement();
+    // NodeValue<Rd> xr = this->_getDisplacement();
 
-    parallel_for(
-      m_given_mesh.numberOfNodes(),
-      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+    // parallel_for(
+    //   m_given_mesh.numberOfNodes(),
+    //   PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
+    //   });
 
-    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    // return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    return nullptr;
   }
 
   std::shared_ptr<const IMesh>
   getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list) const
   {
-    NodeValue<const Rd> given_xr = m_given_mesh.xr();
-
-    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
-
-    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
-    is_displaced.fill(false);
-
-    for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
-      MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
-      const auto cell_list              = cell_zone.cellList();
-      CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
-      is_zone_cell.fill(false);
-      parallel_for(
-        cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
-      parallel_for(
-        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-          auto node_cell_list = node_to_cell_matrix[node_id];
-          bool displace       = true;
-          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
-            const CellId cell_id = node_cell_list[i_node_cell];
-            displace &= is_zone_cell[cell_id];
-          }
-          if (displace) {
-            is_displaced[node_id] = true;
-          }
-        });
-    }
-    synchronize(is_displaced);
-    NodeValue<Rd> xr = this->_getDisplacement();
-
-    parallel_for(
-      m_given_mesh.numberOfNodes(),
-      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
-
-    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    // NodeValue<const Rd> given_xr = m_given_mesh.xr();
+
+    // auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+
+    // NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    // is_displaced.fill(false);
+
+    // for (size_t i_zone = 0; i_zone < zone_descriptor_list.size(); ++i_zone) {
+    //   MeshCellZone<Dimension> cell_zone = getMeshCellZone(m_given_mesh, *zone_descriptor_list[i_zone]);
+    //   const auto cell_list              = cell_zone.cellList();
+    //   CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
+    //   is_zone_cell.fill(false);
+    //   parallel_for(
+    //     cell_list.size(), PUGS_LAMBDA(const size_t i_cell) { is_zone_cell[cell_list[i_cell]] = true; });
+    //   parallel_for(
+    //     m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+    //       auto node_cell_list = node_to_cell_matrix[node_id];
+    //       bool displace       = true;
+    //       for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+    //         const CellId cell_id = node_cell_list[i_node_cell];
+    //         displace &= is_zone_cell[cell_id];
+    //       }
+    //       if (displace) {
+    //         is_displaced[node_id] = true;
+    //       }
+    //     });
+    // }
+    // synchronize(is_displaced);
+    // NodeValue<Rd> xr = this->_getDisplacement();
+
+    // parallel_for(
+    //   m_given_mesh.numberOfNodes(),
+    //   PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
+    //   });
+
+    // return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    return nullptr;
   }
 
   std::shared_ptr<const IMesh>
   getSmoothedMesh(
     const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
   {
-    NodeValue<const Rd> given_xr = m_given_mesh.xr();
+    // NodeValue<const Rd> given_xr = m_given_mesh.xr();
 
-    auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
+    // auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
 
-    NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
-    is_displaced.fill(false);
+    // NodeValue<bool> is_displaced{m_given_mesh.connectivity()};
+    // is_displaced.fill(false);
 
-    for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
-      auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
+    // for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
+    //   auto is_zone_cell = discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const
+    //   double>>();
 
-      parallel_for(
-        m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
-          auto node_cell_list = node_to_cell_matrix[node_id];
-          bool displace       = true;
-          for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
-            const CellId cell_id = node_cell_list[i_node_cell];
-            displace &= (is_zone_cell[cell_id] != 0);
-          }
-          if (displace) {
-            is_displaced[node_id] = true;
-          }
-        });
-    }
-    synchronize(is_displaced);
-    NodeValue<Rd> xr = this->_getDisplacement();
+    //   parallel_for(
+    //     m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
+    //       auto node_cell_list = node_to_cell_matrix[node_id];
+    //       bool displace       = true;
+    //       for (size_t i_node_cell = 0; i_node_cell < node_cell_list.size(); ++i_node_cell) {
+    //         const CellId cell_id = node_cell_list[i_node_cell];
+    //         displace &= (is_zone_cell[cell_id] != 0);
+    //       }
+    //       if (displace) {
+    //         is_displaced[node_id] = true;
+    //       }
+    //     });
+    // }
+    // synchronize(is_displaced);
+    // NodeValue<Rd> xr = this->_getDisplacement();
+
+    // parallel_for(
+    //   m_given_mesh.numberOfNodes(),
+    //   PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id];
+    //   });
 
-    parallel_for(
-      m_given_mesh.numberOfNodes(),
-      PUGS_LAMBDA(const NodeId node_id) { xr[node_id] = is_displaced[node_id] * xr[node_id] + given_xr[node_id]; });
+    // return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
 
-    return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
+    return nullptr;
   }
 
   MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
@@ -519,7 +493,9 @@ class MeshSmootherEscobarHandler::MeshSmootherEscobar
 
   MeshSmootherEscobar(const MeshType& given_mesh,
                       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
-    : m_given_mesh(given_mesh), m_boundary_condition_list(this->_getBCList(given_mesh, bc_descriptor_list))
+    : m_given_mesh{given_mesh},
+      m_boundary_condition_list{this->_getBCList(given_mesh, bc_descriptor_list)},
+      m_is_fixed_node{this->_getFixedNodes()}
   {}
 
   ~MeshSmootherEscobar() = default;