Skip to content
Snippets Groups Projects
Select Git revision
  • cd11ade9ed089f64caa7cf76b9f7a66f1cf53150
  • develop default protected
  • feature/advection
  • feature/composite-scheme-other-fluxes
  • origin/stage/bouguettaia
  • save_clemence
  • feature/local-dt-fsi
  • feature/variational-hydro
  • feature/gmsh-reader
  • feature/reconstruction
  • feature/kinetic-schemes
  • feature/composite-scheme-sources
  • feature/serraille
  • feature/composite-scheme
  • hyperplastic
  • feature/polynomials
  • feature/gks
  • feature/implicit-solver-o2
  • feature/coupling_module
  • feature/implicit-solver
  • feature/merge-local-dt-fsi
  • v0.5.0 protected
  • v0.4.1 protected
  • v0.4.0 protected
  • v0.3.0 protected
  • v0.2.0 protected
  • v0.1.0 protected
  • Kidder
  • v0.0.4 protected
  • v0.0.3 protected
  • v0.0.2 protected
  • v0 protected
  • v0.0.1 protected
33 results

MeshSmootherEscobar.cpp

Blame
  • MeshSmootherEscobar.cpp 54.23 KiB
    #include <mesh/MeshSmootherEscobar.hpp>
    
    #include <algebra/TinyMatrix.hpp>
    #include <algebra/TinyVector.hpp>
    #include <geometry/TriangleTransformation.hpp>
    #include <language/utils/InterpolateItemValue.hpp>
    #include <mesh/Connectivity.hpp>
    #include <mesh/ItemValueUtils.hpp>
    #include <mesh/ItemValueVariant.hpp>
    #include <mesh/Mesh.hpp>
    #include <mesh/MeshCellZone.hpp>
    #include <mesh/MeshFlatNodeBoundary.hpp>
    #include <mesh/MeshLineNodeBoundary.hpp>
    #include <mesh/MeshNodeBoundary.hpp>
    #include <scheme/AxisBoundaryConditionDescriptor.hpp>
    #include <scheme/DiscreteFunctionUtils.hpp>
    #include <scheme/DiscreteFunctionVariant.hpp>
    #include <scheme/FixedBoundaryConditionDescriptor.hpp>
    #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
    #include <utils/RandomEngine.hpp>
    
    #include <variant>
    
    template <size_t Dimension>
    class MeshSmootherEscobarHandler::MeshSmootherEscobar
    {
     private:
      using Rd               = TinyVector<Dimension>;
      using Rdxd             = TinyMatrix<Dimension>;
      using ConnectivityType = Connectivity<Dimension>;
      using MeshType         = Mesh<ConnectivityType>;
    
      const MeshType& m_given_mesh;
      std::shared_ptr<const IMesh> m_p_given_mesh;
    
      class AxisBoundaryCondition;
      class FixedBoundaryCondition;
      class SymmetryBoundaryCondition;
    
      using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
    
      using BoundaryConditionList = std::vector<BoundaryCondition>;
      BoundaryConditionList m_boundary_condition_list;
    
      NodeValue<const bool> m_is_fixed_node;
      NodeValue<const bool> m_is_smoothed_node;
    
      BoundaryConditionList
      _getBCList(const MeshType& mesh,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
      {
        BoundaryConditionList bc_list;
    
        for (const auto& bc_descriptor : bc_descriptor_list) {
          switch (bc_descriptor->type()) {
          case IBoundaryConditionDescriptor::Type::axis: {
            if constexpr (Dimension == 1) {
              bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
            } else {
              bc_list.emplace_back(
                AxisBoundaryCondition{getMeshLineNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
            }
            break;
          }
          case IBoundaryConditionDescriptor::Type::symmetry: {
            bc_list.emplace_back(
              SymmetryBoundaryCondition{getMeshFlatNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
            break;
          }
          case IBoundaryConditionDescriptor::Type::fixed: {
            bc_list.emplace_back(FixedBoundaryCondition{getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())});
            break;
          }
          default: {
            std::ostringstream error_msg;
            error_msg << *bc_descriptor << " is an invalid boundary condition for mesh smoother";
            throw NormalError(error_msg.str());
          }
          }
        }
    
        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);
        }
    
        if constexpr (Dimension == 3) {
          throw NotImplementedError("treat sliding axis node kind");
        }
    
        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(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
      {
        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 Rd& n = bc.outgoingNormal();
    
                const auto node_list = bc.nodeList();
    
                const ConnectivityType& connectivity = m_given_mesh.connectivity();
    
                auto is_boundary_node = connectivity.isBoundaryNode();
                auto node_is_owned    = connectivity.nodeIsOwned();
    
                if constexpr (Dimension == 2) {
                  const TinyVector<Dimension> t{-n[1], n[0]};
    
                  auto cell_to_node_matrix        = connectivity.cellToNodeMatrix();
                  auto node_to_cell_matrix        = connectivity.nodeToCellMatrix();
                  auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
                  auto node_to_edge_matrix        = connectivity.nodeToEdgeMatrix();
                  auto edge_to_node_matrix        = connectivity.edgeToNodeMatrix();
    
                  constexpr double eps = 1E-14;
    
                  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];
    
                    const double alpha = std::acos(-1) / cell_list.size();
    
                    const double cos_alpha     = std::cos(alpha);
                    const double sin_alpha     = std::sin(alpha);
                    const double inv_sin_alpha = 1. / sin_alpha;
                    const double inv_tan_alpha = 1. / std::tan(alpha);
    
                    const TinyMatrix<Dimension> W{1, cos_alpha,   //
                                                  0, sin_alpha};
    
                    const TinyMatrix<Dimension> inv_W = inverse(W);
    
                    TinyMatrix<Dimension, Dimension> dtheta_S =
                      TinyMatrix<Dimension>{t[0], t[0] * (inv_sin_alpha - inv_tan_alpha),   //
                                            t[1], t[1] * (inv_sin_alpha - inv_tan_alpha)};
    
                    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];
                      const 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 = 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]};
    
                      S_list[i_cell] = A * inv_W;
                    }
    
                    SmallArray<double> sigma_list(S_list.size());
                    for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
                      sigma_list[i_cell] = det(S_list[i_cell]);
                    }
    
                    double sigma_min = min(sigma_list);
    
                    double delta = (sigma_min < eps)
                                     ? std::max(std::sqrt(eps * (eps - sigma_min)), std::sqrt(eps) * std::abs(sigma_min))
                                     : 0;
    
                    auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
    
                    double final_f = 0;
    
                    for (size_t i_iter = 0; i_iter < 200; ++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];
                        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 = 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]};
    
                        S_list[i_cell] = A * inv_W;
                      }
    
                      SmallArray<double> sigma_list(S_list.size());
                      for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
                        sigma_list[i_cell] = det(S_list[i_cell]);
                      }
    
                      {
                        double current_sigma_min = min(sigma_list);
                        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;
                        }
                      }
    
                      double f         = 0;
                      double dtheta_f  = 0;
                      double d2theta_f = 0;
    
                      for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
                        const double sigma            = sigma_list[i_cell];
                        const TinyMatrix<Dimension> S = S_list[i_cell];
    
                        const TinyMatrix<Dimension> Sigma = [&S] {
                          TinyMatrix<Dimension> transposed_comS;
                          for (size_t i = 0; i < Dimension; ++i) {
                            for (size_t j = 0; j < Dimension; ++j) {
                              transposed_comS(i, j) = cofactor(S, j, i);
                            }
                          }
                          return transposed_comS;
                        }();
    
                        const double S_norm      = frobenius(S);
                        const double Sigma_norm  = frobenius(Sigma);
                        const double S_norm2     = S_norm * S_norm;
                        const double Sigma_norm2 = Sigma_norm * Sigma_norm;
    
                        const double h = sigma + std::sqrt(sigma * sigma + 4 * delta * delta);
    
                        const double f_jr = S_norm * Sigma_norm / h;
    
                        const double dtheta_sigma = trace(Sigma * dtheta_S);
    
                        TinyMatrix<Dimension, Dimension> dtheta_Sigma =
                          TinyMatrix<Dimension>{t[1] * (inv_sin_alpha - inv_tan_alpha),
                                                -t[0] * (inv_sin_alpha - inv_tan_alpha),   //
                                                -t[1], t[0]};
    
                        const double g = trace(transpose(S) * dtheta_S) / S_norm2                 //
                                         + trace(transpose(Sigma) * dtheta_Sigma) / Sigma_norm2   //
                                         - dtheta_sigma / (h - sigma);
    
                        const double dtheta_f_jr = f_jr * g;
                        const double dtheta2_f_jr =
                          (trace(transpose(dtheta_S) * dtheta_S) / S_norm2                                               //
                           - 2 * trace(transpose(S) * dtheta_S) * trace(transpose(S) * dtheta_S) / (S_norm2 * S_norm2)   //
                                                                                                                         //
                           + trace(transpose(dtheta_Sigma) * dtheta_Sigma) / Sigma_norm2   // + 0
                           - 2 * trace(transpose(Sigma) * dtheta_Sigma) * trace(transpose(Sigma) * dtheta_Sigma) /
                               (Sigma_norm2 * Sigma_norm2)   //
                           //
                           - 2 * trace(dtheta_Sigma * dtheta_S) / (h - sigma)                     //
                           + 2 * sigma / (std::pow(h - sigma, 3)) * dtheta_sigma * dtheta_sigma   //
                           + g * g) *
                          f_jr;
    
                        f += f_jr;
                        dtheta_f += dtheta_f_jr;
                        d2theta_f += dtheta2_f_jr;
                      }
                      if (std::abs(dtheta_f) < 1E-12) {
                        break;
                      }
    
                      const double max_edge_lenght = [=] {
                        auto node_edge_list = node_to_edge_matrix[node_id];
                        double max_lenght   = 0;
                        for (size_t i_edge = 0; i_edge < node_edge_list.size(); ++i_edge) {
                          const EdgeId edge_id            = node_edge_list[i_edge];
                          auto edge_node_list             = edge_to_node_matrix[edge_id];
                          const TinyVector<Dimension>& x0 = old_xr[edge_node_list[0]];
                          const TinyVector<Dimension>& x1 = old_xr[edge_node_list[1]];
                          const double lenght             = l2Norm(x0 - x1);
                          if (lenght > max_lenght) {
                            max_lenght = lenght;
                          }
                        }
                        return max_lenght;
                      }();
    
                      const double delta_x_candidate = dtheta_f / std::abs(d2theta_f);
                      const double delta_x_norm      = std::abs(delta_x_candidate);
    
                      if (delta_x_norm > 0.3 * max_edge_lenght) {
                        x += (0.3 * max_edge_lenght / delta_x_norm) * delta_x_candidate * t;
                      } else {
                        x += delta_x_candidate * t;
                      }
    
                      if (delta_x_norm < 1E-2 * max_edge_lenght) {
                        break;
                      }
    
                      final_f = f;
                    }
                    return final_f;
                  };
    
                  parallel_for(
                    node_list.size(), PUGS_LAMBDA(const size_t i_node) {
                      const NodeId node_id = node_list[i_node];
                      if (m_is_smoothed_node[node_id] and node_is_owned[node_id]) {
                        smooth(node_id, new_xr[node_id]);
                      }
                    });
    
                } else {
                  throw NotImplementedError("Dimension != 2");
                }
    
              } else if constexpr (std::is_same_v<BCType, AxisBoundaryCondition>) {
                if constexpr (Dimension > 1) {
                  throw NotImplementedError("Escobar: axis boundary conditions");
                } else {
                  throw UnexpectedError("AxisBoundaryCondition make no sense in dimension 1");
                }
    
              } else if constexpr (std::is_same_v<BCType, FixedBoundaryCondition>) {
                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];
                    new_xr[node_id]      = old_xr[node_id];
                  });
    
              } else {
                throw UnexpectedError("invalid boundary condition type");
              }
            },
            boundary_condition);
        }
      }
    
      CellValue<const int>
      _getNonConvexCellTag(const NodeValue<const Rd>& xr) const
      {
        CellValue<int> non_convex_cell_tag{m_given_mesh.connectivity()};
    
        auto cell_to_node_matrix = m_given_mesh.connectivity().cellToNodeMatrix();
    
        if constexpr (Dimension == 2) {
          parallel_for(
            m_given_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
              auto cell_node_list = cell_to_node_matrix[cell_id];
              bool is_convex      = true;
              for (size_t i_node = 0; i_node < cell_node_list.size(); ++i_node) {
                const TriangleTransformation<Dimension> T{xr[cell_node_list[i_node]],
                                                          xr[cell_node_list[(i_node + 1) % cell_node_list.size()]],
                                                          xr[cell_node_list[(i_node + 2) % cell_node_list.size()]]};
                if (T.jacobianDeterminant() <= 0) {
                  is_convex = false;
                  break;
                }
              }
              non_convex_cell_tag[cell_id] = (is_convex == false);
            });
        } else {
          throw NotImplementedError("non convex cell calculation in dimension != 2");
        }
    
        return non_convex_cell_tag;
      }
    
      NodeValue<const bool>
      _getSmoothedNode(CellValue<const int>& non_convex_cell_tag, const NodeValue<const bool>& m_is_fixed_node) const
      {
        auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
        NodeValue<bool> is_smoothed{m_given_mesh.connectivity()};
        is_smoothed.fill(false);
        parallel_for(
          m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
            if (not m_is_fixed_node[node_id]) {
              auto node_cell_list = node_to_cell_matrix[node_id];
              for (size_t i_cell = 0; i_cell < node_cell_list.size(); ++i_cell) {
                const CellId cell_id = node_cell_list[i_cell];
                if (non_convex_cell_tag[cell_id]) {
                  is_smoothed[node_id] = true;
                  break;
                }
              }
            }
          });
    
        return is_smoothed;
      }
    
     public:
      std::shared_ptr<const ItemValueVariant>
      getQuality() const
      {
        return std::make_shared<ItemValueVariant>(m_given_mesh.xr());
      }
    
      void
      _getNewPositions(const NodeValue<const Rd>& old_xr, NodeValue<Rd>& new_xr) const
      {
        constexpr bool multi_connection = true;
    
        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_to_edge_matrix        = connectivity.nodeToEdgeMatrix();
          auto edge_to_node_matrix        = connectivity.edgeToNodeMatrix();
          auto node_number_in_their_cells = connectivity.nodeLocalNumbersInTheirCells();
    
          NodeValue<double> quality{connectivity};
    
          constexpr double eps = 1E-14;
          quality.fill(2);
    
          auto smooth = [=](const NodeId node_id, TinyVector<Dimension>& x) -> double {
            const bool print = (node_id == 4);
    
            auto cell_list           = node_to_cell_matrix[node_id];
            auto node_number_in_cell = node_number_in_their_cells[node_id];
    
            const double alpha         = 2 * std::acos(-1) / cell_list.size();
            const double sin_alpha     = std::sin(alpha);
            const double cos_alpha     = std::cos(alpha);
            const double inv_sin_alpha = 1. / sin_alpha;
            const double inv_tan_alpha = 1. / std::tan(alpha);
    
            const TinyMatrix<Dimension> W{1, cos_alpha,   //
                                          0, sin_alpha};
    
            const TinyMatrix<Dimension> inv_W = inverse(W);
    
            std::array<TinyMatrix<Dimension>, Dimension>    //
              S_gradient = {TinyMatrix<Dimension>{-1, -1,   //
                                                  +0, +0},
                            TinyMatrix<Dimension>{+inv_tan_alpha, +inv_tan_alpha,   //
                                                  -inv_sin_alpha, -inv_sin_alpha}};
    
            SmallArray<TinyMatrix<Dimension>> S_list(cell_list.size());
    
            double sigma_min = std::numeric_limits<double>::max();
            double delta     = 0;
    
            auto frobenius = [](const TinyMatrix<Dimension>& M) { return std::sqrt(trace(transpose(M) * M)); };
    
            double final_f = 0;
    
            std::cout.precision(15);
    
            for (size_t i_iter = 0; i_iter < 1; ++i_iter) {
              if (print) {
                std::cout << "------------------------------\n";
                std::cout << "iter = " << i_iter << '\n';
                std::cout << "x    = " << x << '\n';
              }
    
              double f                         = 0;
              TinyVector<Dimension> f_gradient = zero;
              TinyMatrix<Dimension> f_hessian  = zero;
    
              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];
                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 = 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]};
    
                S_list[i_cell] = A * inv_W;
    
                if (print) {
                  std::cout << "--- S --- cell = " << i_cell << " (" << node_to_cell_matrix[node_id][i_cell] << ")\n";
                  std::cout << "a = " << a << '\n';
                  std::cout << "b = " << b << '\n';
                  std::cout << "x = " << x << '\n';
                  std::cout << "A = " << A << '\n';
                  std::cout << "S = " << S_list[i_cell] << '\n';
                }
              }
    
              SmallArray<double> sigma_list(S_list.size());
              for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
                sigma_list[i_cell] = det(S_list[i_cell]);
              }
    
              {
                double current_sigma_min = min(sigma_list);
    
                if constexpr (multi_connection) {
                  size_t nb_medium_range_simplices = 0;
                  for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
                    const CellId cell_id = cell_list[i_cell];
                    nb_medium_range_simplices += 2 * (cell_to_node_matrix[cell_id].size() - 3);
                  }
    
                  if (nb_medium_range_simplices > 0) {
                    SmallArray<double> medium_range_sigma_list(nb_medium_range_simplices);
    
                    for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
                      const CellId cell_id       = cell_list[i_cell];
                      const size_t i_cell_node   = node_number_in_cell[i_cell];
                      auto cell_node_list        = cell_to_node_matrix[cell_list[i_cell]];
                      const size_t cell_nb_nodes = cell_node_list.size();
    
                      const double mr_alpha     = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
                      const double mr_sin_alpha = std::sin(mr_alpha);
                      const double mr_cos_alpha = std::cos(mr_alpha);
    
                      const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha,   //
                                                       0, mr_sin_alpha};
    
                      const TinyMatrix mr_inv_W = inverse(mr_W);
    
                      if (print) {
                        std::cout << "MC === i_cell = " << i_cell << '\n';
                      }
    
                      {   // first simplex
                        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 + 2) % cell_nb_nodes]];
    
                        const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0],   //
                                                         a[1] - x[1], b[1] - x[1]};
    
                        const TinyMatrix mr_S = mr_A * mr_inv_W;
                        const double mr_sigma = det(mr_S);
                        if (mr_sigma < current_sigma_min) {
                          current_sigma_min = mr_sigma;
                        }
    
                        if (print) {
                          std::cout << "MC triangle 1 ===\n";
                          std::cout << "MC a = " << a << '\n';
                          std::cout << "MC b = " << b << '\n';
                          std::cout << "MC x = " << x << '\n';
                        }
                      }
    
                      {   // second simplex
                        const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
                        const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
    
                        const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0],   //
                                                         a[1] - x[1], b[1] - x[1]};
    
                        const TinyMatrix mr_S = mr_A * mr_inv_W;
                        const double mr_sigma = det(mr_S);
                        if (mr_sigma < current_sigma_min) {
                          current_sigma_min = mr_sigma;
                        }
    
                        if (print) {
                          std::cout << "MC triangle 2 ===\n";
                          std::cout << "MC a = " << a << '\n';
                          std::cout << "MC b = " << b << '\n';
                          std::cout << "MC x = " << x << '\n';
                        }
                      }
                    }
                  }
                }
    
                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;
                }
    
                if (print) {
                  std::cout << " ==== \n";
                  std::cout << "delta = " << delta << '\n';
                }
              }
    
              auto escobarSSigma = [&frobenius, delta,   //
                                                         //
                                    &print, &node_to_cell_matrix,
                                    &node_id](const double sigma, const TinyMatrix<Dimension>& S,
                                              const TinyVector<Dimension>& sigma_gradient,
                                              const std::array<TinyMatrix<Dimension>, Dimension>& S_gradient,
                                              const std::array<TinyMatrix<Dimension>, Dimension>& Sigma_gradient, double& f,
                                              TinyVector<Dimension>& f_gradient, TinyMatrix<Dimension>& f_hessian) {
                const TinyMatrix<Dimension> Sigma = [&S] {
                  TinyMatrix<Dimension> transposed_comS;
                  for (size_t i = 0; i < Dimension; ++i) {
                    for (size_t j = 0; j < Dimension; ++j) {
                      transposed_comS(i, j) = cofactor(S, j, i);
                    }
                  }
                  return transposed_comS;
                }();
    
                const double S_norm      = frobenius(S);
                const double Sigma_norm  = frobenius(Sigma);
                const double S_norm2     = S_norm * S_norm;
                const double Sigma_norm2 = Sigma_norm * Sigma_norm;
    
                double h = sigma + std::sqrt(sigma * sigma + 4 * delta * delta);
    
                const double f_jr = S_norm * Sigma_norm / h;
    
                TinyVector<Dimension> g{trace(transpose(S) * S_gradient[0]) / S_norm2                   //
                                          + trace(transpose(Sigma) * Sigma_gradient[0]) / Sigma_norm2   //
                                          - trace(Sigma * S_gradient[0]) / (h - sigma),
                                        //
                                        trace(transpose(S) * S_gradient[1]) / S_norm2                   //
                                          + trace(transpose(Sigma) * Sigma_gradient[1]) / Sigma_norm2   //
                                          - trace(Sigma * S_gradient[1]) / (h - sigma)};
    
                const TinyVector<Dimension> f_jr_gradient = f_jr * g;
    
                TinyMatrix<Dimension> f_jr_hessian = zero;
                for (size_t i = 0; i < Dimension; ++i) {
                  for (size_t j = 0; j < Dimension; ++j) {
                    f_jr_hessian(i, j) =                                           //
                      (trace(transpose(S_gradient[j]) * S_gradient[i]) / S_norm2   //
                       - 2 * trace(transpose(S) * S_gradient[j]) * trace(transpose(S) * S_gradient[i]) /
                           (S_norm2 * S_norm2)                                                   //
                                                                                                 //
                       + 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)   //
    
                       - trace(Sigma_gradient[j] * S_gradient[i]) / (h - sigma)                     //
                       + sigma / (std::pow(h - sigma, 3)) * sigma_gradient[i] * sigma_gradient[j]   //
                       + g[j] * g[i]) *
                      f_jr;
                  }
                }
    
                if (print) {
                  std::cout << "S     = " << S << '\n';
                  std::cout << "Sigma = " << Sigma << '\n';
    
                  std::cout << "normS     = " << S_norm << '\n';
                  std::cout << "normSigma = " << Sigma_norm << '\n';
    
                  std::cout << "sigma = " << sigma << '\n';
                  std::cout << "f_jr  = " << f_jr << '\n';
                  std::cout << "g     = " << g << '\n';
                  std::cout << "grad_f_jr = " << f_jr_gradient << '\n';
                  std::cout << "hess_f_jr = " << f_jr_hessian << '\n';
                }
    
                f += f_jr;
                f_gradient += f_jr_gradient;
                f_hessian += f_jr_hessian;
              };
    
              auto escobar = escobarSSigma;
    
              for (size_t i_cell = 0; i_cell < S_list.size(); ++i_cell) {
                if (print) {
                  std::cout << "------ cell = " << i_cell << '\n';
                }
    
                const double sigma            = sigma_list[i_cell];
                const TinyMatrix<Dimension> S = S_list[i_cell];
    
                const TinyMatrix<Dimension> Sigma = [&S] {
                  TinyMatrix<Dimension> transposed_comS;
                  for (size_t i = 0; i < Dimension; ++i) {
                    for (size_t j = 0; j < Dimension; ++j) {
                      transposed_comS(i, j) = cofactor(S, j, i);
                    }
                  }
                  return transposed_comS;
                }();
    
                TinyVector<Dimension> sigma_gradient{trace(Sigma * S_gradient[0]),   //
                                                     trace(Sigma * S_gradient[1])};
    
                const std::array<TinyMatrix<Dimension>, Dimension>   //
                  Sigma_gradient{TinyMatrix<Dimension>{+0, +1,       //
                                                       +0, -1},
                                 TinyMatrix<Dimension>{-inv_sin_alpha, -inv_tan_alpha,   //
                                                       +inv_sin_alpha, +inv_tan_alpha}};
    
                escobar(sigma, S, sigma_gradient, S_gradient, Sigma_gradient, f, f_gradient, f_hessian);
              }
    
              if constexpr (multi_connection) {
                for (size_t i_cell = 0; i_cell < cell_list.size(); ++i_cell) {
                  const CellId cell_id       = cell_list[i_cell];
                  auto cell_node_list        = cell_to_node_matrix[cell_list[i_cell]];
                  const size_t cell_nb_nodes = cell_node_list.size();
    
                  if (print) {
                    std::cout << "=== MC cell " << i_cell << " ===\n";
                  }
    
                  if (cell_nb_nodes <= 3) {
                    continue;
                  }
    
                  const size_t i_cell_node = node_number_in_cell[i_cell];
    
                  const double mr_alpha = 2 * std::acos(-1) / ((cell_nb_nodes - 2) * (cell_list.size()));
                  if (print) {
                    std::cout << "=== MC mr_alpha " << mr_alpha << " ===\n";
                  }
                  const double mr_sin_alpha     = std::sin(mr_alpha);
                  const double mr_cos_alpha     = std::cos(mr_alpha);
                  const double mr_inv_sin_alpha = 1. / mr_sin_alpha;
                  const double mr_inv_tan_alpha = 1. / std::tan(mr_alpha);
    
                  const TinyMatrix<Dimension> mr_W{1, mr_cos_alpha,   //
                                                   0, mr_sin_alpha};
    
                  const TinyMatrix mr_inv_W = inverse(mr_W);
    
                  std::array<TinyMatrix<Dimension>, Dimension>                                      //
                    mr_S_gradient = {TinyMatrix<Dimension>{-1, (mr_cos_alpha - 1) / mr_sin_alpha,   //
                                                           +0, +0},
                                     TinyMatrix<Dimension>{+0, +0,   //
                                                           -1, (mr_cos_alpha - 1) / mr_sin_alpha}};
    
                  {   // first simplex
                    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 + 2) % cell_nb_nodes]];
    
                    const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0],   //
                                                     a[1] - x[1], b[1] - x[1]};
    
                    const TinyMatrix mr_S = mr_A * mr_inv_W;
                    const double mr_sigma = det(mr_S);
    
                    const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
                      TinyMatrix<Dimension> transposed_comS;
                      for (size_t i = 0; i < Dimension; ++i) {
                        for (size_t j = 0; j < Dimension; ++j) {
                          transposed_comS(i, j) = cofactor(mr_S, j, i);
                        }
                      }
                      return transposed_comS;
                    }();
    
                    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 - 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 1 ===\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);
                  }
    
                  {   // second simplex
                    const TinyVector a = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 2) % cell_nb_nodes]];
                    const TinyVector b = old_xr[cell_node_list[(i_cell_node + cell_nb_nodes - 1) % cell_nb_nodes]];
    
                    const TinyMatrix<Dimension> mr_A{a[0] - x[0], b[0] - x[0],   //
                                                     a[1] - x[1], b[1] - x[1]};
    
                    const TinyMatrix mr_S = mr_A * mr_inv_W;
                    const double mr_sigma = det(mr_S);
    
                    const TinyMatrix<Dimension> mr_Sigma = [&mr_S] {
                      TinyMatrix<Dimension> transposed_comS;
                      for (size_t i = 0; i < Dimension; ++i) {
                        for (size_t j = 0; j < Dimension; ++j) {
                          transposed_comS(i, j) = cofactor(mr_S, j, i);
                        }
                      }
                      return transposed_comS;
                    }();
    
                    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}};
    
                    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 << " mr_S = " << mr_S << '\n';
                      std::cout << " mr_Sigma = " << mr_Sigma << '\n';
                    }
    
                    escobar(mr_sigma, mr_S, mr_sigma_gradient, mr_S_gradient, mr_Sigma_gradient, f, f_gradient, f_hessian);
                  }
                }
              }
    
              const double max_edge_lenght = [=] {
                auto node_edge_list = node_to_edge_matrix[node_id];
                double max_lenght   = 0;
                for (size_t i_edge = 0; i_edge < node_edge_list.size(); ++i_edge) {
                  const EdgeId edge_id            = node_edge_list[i_edge];
                  auto edge_node_list             = edge_to_node_matrix[edge_id];
                  const TinyVector<Dimension>& x0 = old_xr[edge_node_list[0]];
                  const TinyVector<Dimension>& x1 = old_xr[edge_node_list[1]];
                  const double lenght             = l2Norm(x0 - x1);
                  if (lenght > max_lenght) {
                    max_lenght = lenght;
                  }
                }
                return max_lenght;
              }();
    
              const TinyVector delta_x_candidate = -inverse(f_hessian) * f_gradient;
    
              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 << "f_gradient = " << f_gradient << '\n';
                std::cout << "f_hessian  = " << f_hessian << '\n';
                std::cout << "inv(f_hes) = " << inverse(f_hessian) << '\n';
                std::cout << "factor = " << factor << '\n';
                std::cout << "delta_x= " << delta_x_candidate << '\n';
              }
              if (l2Norm(f_gradient) < 1E-12) {
                break;
              }
    
              if (delta_x_norm < 1E-8 * max_edge_lenght) {
                break;
              }
    
              final_f = f;
            }
    
            return final_f;
          };
    
          // parallel_for(
          //   connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId node_id) {
          //     if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
          //       quality[node_id] = smooth(node_id, new_xr[node_id]);
          //     }
          //   });
    
          size_t count = 0;
          for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
            if (node_is_owned[node_id] and m_is_smoothed_node[node_id] and not is_boundary_node[node_id]) {
              quality[node_id] = smooth(node_id, new_xr[node_id]);
              count++;
            }
          }
          std::cout << "nb smoothed = " << count << '\n';
        } else {
          throw NotImplementedError("Dimension != 2");
        }
        this->_applyBC(old_xr, new_xr);
      }
    
      std::shared_ptr<const IMesh>
      _smoothedMesh()
      {
        NodeValue<const Rd> xr = m_given_mesh.xr();
    
        size_t nb_non_convex_cells = 0;
        size_t i_convexify         = 1;
        do {
          std::cout << "- convexfication step " << i_convexify++ << '\n';
          NodeValue<Rd> new_xr = copy(xr);
          this->_getNewPositions(xr, new_xr);
          xr                       = new_xr;
          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) {
            m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
          }
        } while ((nb_non_convex_cells > 0) and (i_convexify < 10));
    
        return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh()
      {
        NodeValue<bool> is_smoothed_node(m_given_mesh.connectivity());
        parallel_for(
          m_given_mesh.numberOfNodes(),
          PUGS_LAMBDA(const NodeId node_id) { is_smoothed_node[node_id] = not m_is_fixed_node[node_id]; });
    
        m_is_smoothed_node = is_smoothed_node;
    
        return _smoothedMesh();
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh(const FunctionSymbolId& function_symbol_id)
      {
        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<bool> is_fixed_node(m_given_mesh.connectivity());
        parallel_for(
          m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
            is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
          });
    
        synchronize(is_fixed_node);
    
        m_is_fixed_node = is_fixed_node;
    
        return _smoothedMesh();
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh(const std::vector<std::shared_ptr<const IZoneDescriptor>>& zone_descriptor_list)
      {
        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<bool> is_fixed_node(m_given_mesh.connectivity());
        parallel_for(
          m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
            is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
          });
    
        synchronize(is_fixed_node);
    
        m_is_fixed_node = is_fixed_node;
    
        return _smoothedMesh();
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh(const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list)
      {
        NodeValue<const Rd> given_xr = m_given_mesh.xr();
    
        auto node_to_cell_matrix = m_given_mesh.connectivity().nodeToCellMatrix();
    
        CellValue<bool> is_zone_cell{m_given_mesh.connectivity()};
        is_zone_cell.fill(false);
    
        for (size_t i_zone = 0; i_zone < discrete_function_variant_list.size(); ++i_zone) {
          auto characteristic_function =
            discrete_function_variant_list[i_zone]->get<DiscreteFunctionP0<Dimension, const double>>();
          parallel_for(
            m_given_mesh.numberOfCells(), PUGS_LAMBDA(const CellId cell_id) {
              if (characteristic_function[cell_id] != 0) {
                is_zone_cell[cell_id] = true;
              }
            });
        }
    
        CellValue<const int> non_convex_cell_tag = this->_getNonConvexCellTag(given_xr);
    
        const bool has_non_convex_cell = max(non_convex_cell_tag) > 0;
    
        if (not has_non_convex_cell) {
          return m_p_given_mesh;
        }
    
        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 characteristic_function =
            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 &= (characteristic_function[cell_id] != 0);
              }
              if (displace) {
                is_displaced[node_id] = true;
              }
            });
        }
        synchronize(is_displaced);
    
        NodeValue<bool> is_fixed_node(m_given_mesh.connectivity());
        parallel_for(
          m_given_mesh.numberOfNodes(), PUGS_LAMBDA(const NodeId node_id) {
            is_fixed_node[node_id] = m_is_fixed_node[node_id] or not is_displaced[node_id];
          });
    
        synchronize(is_fixed_node);
    
        m_is_fixed_node = is_fixed_node;
    
        m_is_smoothed_node = _getSmoothedNode(non_convex_cell_tag, m_is_fixed_node);
    
        return _smoothedMesh();
      }
    
      MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
      MeshSmootherEscobar(MeshSmootherEscobar&&)      = delete;
    
      MeshSmootherEscobar(const MeshType& given_mesh,
                          const std::shared_ptr<const IMesh>& p_mesh,
                          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
        : m_given_mesh{given_mesh},
          m_p_given_mesh{p_mesh},
          m_boundary_condition_list{this->_getBCList(given_mesh, bc_descriptor_list)},
          m_is_fixed_node{this->_getFixedNodes()}
      {}
    
      ~MeshSmootherEscobar() = default;
    };
    
    template <size_t Dimension>
    class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::AxisBoundaryCondition
    {
     public:
      using Rd = TinyVector<Dimension, double>;
    
     private:
      const MeshLineNodeBoundary<Dimension> m_mesh_line_node_boundary;
    
     public:
      const Rd&
      direction() const
      {
        return m_mesh_line_node_boundary.direction();
      }
    
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_line_node_boundary.nodeList();
      }
    
      AxisBoundaryCondition(MeshLineNodeBoundary<Dimension>&& mesh_line_node_boundary)
        : m_mesh_line_node_boundary(mesh_line_node_boundary)
      {
        ;
      }
    
      ~AxisBoundaryCondition() = default;
    };
    
    template <>
    class MeshSmootherEscobarHandler::MeshSmootherEscobar<1>::AxisBoundaryCondition
    {
     public:
      AxisBoundaryCondition()  = default;
      ~AxisBoundaryCondition() = default;
    };
    
    template <size_t Dimension>
    class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::FixedBoundaryCondition
    {
     private:
      const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
    
     public:
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_node_boundary.nodeList();
      }
    
      FixedBoundaryCondition(MeshNodeBoundary<Dimension>&& mesh_node_boundary) : m_mesh_node_boundary{mesh_node_boundary} {}
    
      ~FixedBoundaryCondition() = default;
    };
    
    template <size_t Dimension>
    class MeshSmootherEscobarHandler::MeshSmootherEscobar<Dimension>::SymmetryBoundaryCondition
    {
     public:
      using Rd = TinyVector<Dimension, double>;
    
     private:
      const MeshFlatNodeBoundary<Dimension> m_mesh_flat_node_boundary;
    
     public:
      const Rd&
      outgoingNormal() const
      {
        return m_mesh_flat_node_boundary.outgoingNormal();
      }
    
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_flat_node_boundary.nodeList();
      }
    
      SymmetryBoundaryCondition(MeshFlatNodeBoundary<Dimension>&& mesh_flat_node_boundary)
        : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
      {
        ;
      }
    
      ~SymmetryBoundaryCondition() = default;
    };
    
    std::shared_ptr<const ItemValueVariant>
    MeshSmootherEscobarHandler::getQuality(
      const std::shared_ptr<const IMesh>& mesh,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
    {
      switch (mesh->dimension()) {
      case 1: {
        constexpr size_t Dimension = 1;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getQuality();
      }
      case 2: {
        constexpr size_t Dimension = 2;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getQuality();
      }
      case 3: {
        constexpr size_t Dimension = 3;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getQuality();
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }
    
    std::shared_ptr<const IMesh>
    MeshSmootherEscobarHandler::getSmoothedMesh(
      const std::shared_ptr<const IMesh>& mesh,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
    {
      switch (mesh->dimension()) {
      case 1: {
        constexpr size_t Dimension = 1;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh();
      }
      case 2: {
        constexpr size_t Dimension = 2;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh();
      }
      case 3: {
        constexpr size_t Dimension = 3;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh();
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }
    
    std::shared_ptr<const IMesh>
    MeshSmootherEscobarHandler::getSmoothedMesh(
      const std::shared_ptr<const IMesh>& mesh,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
      const FunctionSymbolId& function_symbol_id) const
    {
      switch (mesh->dimension()) {
      case 1: {
        constexpr size_t Dimension = 1;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(function_symbol_id);
      }
      case 2: {
        constexpr size_t Dimension = 2;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(function_symbol_id);
      }
      case 3: {
        constexpr size_t Dimension = 3;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(function_symbol_id);
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }
    
    std::shared_ptr<const IMesh>
    MeshSmootherEscobarHandler::getSmoothedMesh(
      const std::shared_ptr<const IMesh>& mesh,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
      const std::vector<std::shared_ptr<const IZoneDescriptor>>& smoothing_zone_list) const
    {
      switch (mesh->dimension()) {
      case 1: {
        constexpr size_t Dimension = 1;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(smoothing_zone_list);
      }
      case 2: {
        constexpr size_t Dimension = 2;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(smoothing_zone_list);
      }
      case 3: {
        constexpr size_t Dimension = 3;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(smoothing_zone_list);
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }
    
    std::shared_ptr<const IMesh>
    MeshSmootherEscobarHandler::getSmoothedMesh(
      const std::shared_ptr<const IMesh>& mesh,
      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
      const std::vector<std::shared_ptr<const DiscreteFunctionVariant>>& discrete_function_variant_list) const
    {
      if (not hasSameMesh(discrete_function_variant_list)) {
        throw NormalError("discrete functions are not defined on the same mesh");
      }
    
      std::shared_ptr<const IMesh> common_mesh = getCommonMesh(discrete_function_variant_list);
    
      if (common_mesh != mesh) {
        throw NormalError("discrete functions are not defined on the smoothed mesh");
      }
    
      switch (mesh->dimension()) {
      case 1: {
        constexpr size_t Dimension = 1;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(discrete_function_variant_list);
      }
      case 2: {
        constexpr size_t Dimension = 2;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(discrete_function_variant_list);
      }
      case 3: {
        constexpr size_t Dimension = 3;
        using MeshType             = Mesh<Connectivity<Dimension>>;
        MeshSmootherEscobar smoother(dynamic_cast<const MeshType&>(*mesh), mesh, bc_descriptor_list);
        return smoother.getSmoothedMesh(discrete_function_variant_list);
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }