Skip to content
Snippets Groups Projects
Select Git revision
  • c9d357071e3a9be7d9e5d3bcd915673687bb23b5
  • develop default protected
  • feature/variational-hydro
  • origin/stage/bouguettaia
  • feature/gmsh-reader
  • feature/reconstruction
  • save_clemence
  • feature/kinetic-schemes
  • feature/local-dt-fsi
  • feature/composite-scheme-sources
  • feature/composite-scheme-other-fluxes
  • 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
  • master protected
  • 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

ASTBuilder.cpp

Blame
  • MeshSmootherEscobar.cpp 28.93 KiB
    #include <mesh/MeshSmootherEscobar.hpp>
    
    #include <algebra/TinyMatrix.hpp>
    #include <algebra/TinyVector.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;
    
      class AxisBoundaryCondition;
      class FixedBoundaryCondition;
      class SymmetryBoundaryCondition;
    
      using BoundaryCondition = std::variant<AxisBoundaryCondition, FixedBoundaryCondition, SymmetryBoundaryCondition>;
    
      using BoundaryConditionList = std::vector<BoundaryCondition>;
      BoundaryConditionList m_boundary_condition_list;
    
      BoundaryConditionList
      _getBCList(const MeshType& mesh,
                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list)
      {
        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;
      }
    
      void
      _applyBC(NodeValue<Rd>& shift) 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 Rdxd I   = identity;
                const Rdxd nxn = tensorProduct(n, n);
                const Rdxd P   = I - nxn;
    
                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] = P * shift[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];
                    });
                } 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];
                    shift[node_id]       = zero;
                  });
    
              } else {
                throw UnexpectedError("invalid boundary condition type");
              }
            },
            boundary_condition);
        }
      }
    
      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();
    
          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 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 TinyMatrix<Dimension> W{1, std::cos(alpha),   //
                                          0, std::sin(alpha)};
    
            const TinyMatrix<Dimension> inv_W = inverse(W);
    
            std::array<TinyMatrix<Dimension>, Dimension> S_gradient =
              {TinyMatrix<Dimension>{-1, -1. / std::sin(alpha) + 1. / std::tan(alpha),   //
                                     +0, +0},                                            //
               TinyMatrix<Dimension>{+0, +0,                                             //
                                     -1, -1. / std::sin(alpha) + 1. / std::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];
              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 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]);
            }
    
            const double sigma_min = min(sigma_list);
    
            const 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)); };
    
            // 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) {
              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 = 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 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 f                         = 0;
              TinyVector<Dimension> f_gradient = zero;
              TinyMatrix<Dimension> f_hessian  = zero;
    
              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 = sigma * inverse(S);
    
                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;
    
                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;
    
                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)   //
                       //
                       - 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]) *
                      f_jr;
                  }
                }
    
                f += f_jr;
                f_gradient += f_jr_gradient;
                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;
              }
    
              x -= inverse(f_hessian) * f_gradient;
    
              final_f = f;
            }
            return final_f;
          };
    
          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;
              }
            });
    
          return std::make_shared<ItemValueVariant>(quality);
        } else {
          throw NotImplementedError("Dimension != 2");
        }
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh() const
      {
        NodeValue<const Rd> given_xr = m_given_mesh.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]; });
    
        return std::make_shared<MeshType>(m_given_mesh.shared_connectivity(), xr);
      }
    
      std::shared_ptr<const IMesh>
      getSmoothedMesh(const FunctionSymbolId& function_symbol_id) const
      {
        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<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);
      }
    
      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);
      }
    
      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();
    
        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 < 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) { 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);
      }
    
      MeshSmootherEscobar(const MeshSmootherEscobar&) = delete;
      MeshSmootherEscobar(MeshSmootherEscobar&&)      = delete;
    
      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))
      {}
    
      ~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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), 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), bc_descriptor_list);
        return smoother.getSmoothedMesh(discrete_function_variant_list);
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }