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

IUnaryOperatorProcessorBuilder.hpp

Blame
  • HyperelasticSolver.cpp 30.81 KiB
    #include <scheme/HyperelasticSolver.hpp>
    
    #include <language/utils/InterpolateItemValue.hpp>
    #include <mesh/ItemValueUtils.hpp>
    #include <mesh/ItemValueVariant.hpp>
    #include <mesh/MeshFaceBoundary.hpp>
    #include <mesh/MeshFlatNodeBoundary.hpp>
    #include <mesh/MeshNodeBoundary.hpp>
    #include <mesh/SubItemValuePerItemVariant.hpp>
    #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
    #include <scheme/DiscreteFunctionP0.hpp>
    #include <scheme/DiscreteFunctionUtils.hpp>
    #include <scheme/ExternalBoundaryConditionDescriptor.hpp>
    #include <scheme/FixedBoundaryConditionDescriptor.hpp>
    #include <scheme/IBoundaryConditionDescriptor.hpp>
    #include <scheme/IDiscreteFunction.hpp>
    #include <scheme/IDiscreteFunctionDescriptor.hpp>
    #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
    #include <utils/Socket.hpp>
    
    #include <variant>
    #include <vector>
    
    template <size_t Dimension>
    double
    hyperelastic_dt(const DiscreteFunctionP0<Dimension, double>& c)
    {
      const Mesh<Connectivity<Dimension>>& mesh = dynamic_cast<const Mesh<Connectivity<Dimension>>&>(*c.mesh());
    
      const auto Vj = MeshDataManager::instance().getMeshData(mesh).Vj();
      const auto Sj = MeshDataManager::instance().getMeshData(mesh).sumOverRLjr();
    
      CellValue<double> local_dt{mesh.connectivity()};
      parallel_for(
        mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { local_dt[j] = 2 * Vj[j] / (Sj[j] * c[j]); });
    
      return min(local_dt);
    }
    
    double
    hyperelastic_dt(const std::shared_ptr<const IDiscreteFunction>& c)
    {
      if ((c->descriptor().type() != DiscreteFunctionType::P0) or (c->dataType() != ASTNodeDataType::double_t)) {
        throw NormalError("invalid discrete function type");
      }
    
      std::shared_ptr mesh = c->mesh();
    
      switch (mesh->dimension()) {
      case 1: {
        return hyperelastic_dt(dynamic_cast<const DiscreteFunctionP0<1, double>&>(*c));
      }
      case 2: {
        return hyperelastic_dt(dynamic_cast<const DiscreteFunctionP0<2, double>&>(*c));
      }
      case 3: {
        return hyperelastic_dt(dynamic_cast<const DiscreteFunctionP0<3, double>&>(*c));
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }
    
    template <size_t Dimension>
    class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticSolverHandler::IHyperelasticSolver
    {
     private:
      using Rdxd = TinyMatrix<Dimension>;
      using Rd   = TinyVector<Dimension>;
    
      using MeshType     = Mesh<Connectivity<Dimension>>;
      using MeshDataType = MeshData<Dimension>;
    
      using DiscreteScalarFunction = DiscreteFunctionP0<Dimension, double>;
      using DiscreteVectorFunction = DiscreteFunctionP0<Dimension, Rd>;
      using DiscreteTensorFunction = DiscreteFunctionP0<Dimension, Rdxd>;
    
      class FixedBoundaryCondition;
      class PressureBoundaryCondition;
      class NormalStressBoundaryCondition;
      class SymmetryBoundaryCondition;
      class VelocityBoundaryCondition;
    
      using BoundaryCondition = std::
        variant<FixedBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
    
      using BoundaryConditionList = std::vector<BoundaryCondition>;
    
      NodeValuePerCell<const Rdxd>
      _computeGlaceAjr(const MeshType& mesh, const DiscreteScalarFunction& rhoaL, const DiscreteScalarFunction& rhoaT) const
      {
        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
    
        const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
        const NodeValuePerCell<const Rd> njr = mesh_data.njr();
    
        NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
        const Rdxd I = identity;
        parallel_for(
          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
            const size_t& nb_nodes = Ajr.numberOfSubValues(j);
            const double& rhoaL_j  = rhoaL[j];
            const double& rhoaT_j  = rhoaT[j];
            for (size_t r = 0; r < nb_nodes; ++r) {
              const Rdxd& M = tensorProduct(Cjr(j, r), njr(j, r));
              Ajr(j, r)     = rhoaL_j * M + rhoaT_j * (I - M);
            }
          });
    
        return Ajr;
      }
    
      NodeValuePerCell<const Rdxd>
      _computeEucclhydAjr(const MeshType& mesh,
                          const DiscreteScalarFunction& rhoaL,
                          const DiscreteScalarFunction& rhoaT) const
      {
        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
    
        const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
        const NodeValuePerFace<const Rd> nlr = mesh_data.nlr();
    
        const auto& face_to_node_matrix = mesh.connectivity().faceToNodeMatrix();
        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
        const auto& cell_to_face_matrix = mesh.connectivity().cellToFaceMatrix();
    
        NodeValuePerCell<Rdxd> Ajr{mesh.connectivity()};
    
        parallel_for(
          Ajr.numberOfValues(), PUGS_LAMBDA(size_t jr) { Ajr[jr] = zero; });
        const Rdxd I = identity;
        parallel_for(
          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
            const auto& cell_nodes = cell_to_node_matrix[j];
    
            const auto& cell_faces = cell_to_face_matrix[j];
    
            const double& rho_aL = rhoaL[j];
            const double& rho_aT = rhoaT[j];
    
            for (size_t L = 0; L < cell_faces.size(); ++L) {
              const FaceId& l        = cell_faces[L];
              const auto& face_nodes = face_to_node_matrix[l];
    
              auto local_node_number_in_cell = [&](NodeId node_number) {
                for (size_t i_node = 0; i_node < cell_nodes.size(); ++i_node) {
                  if (node_number == cell_nodes[i_node]) {
                    return i_node;
                  }
                }
                return std::numeric_limits<size_t>::max();
              };
    
              for (size_t rl = 0; rl < face_nodes.size(); ++rl) {
                const size_t R = local_node_number_in_cell(face_nodes[rl]);
                const Rdxd& M  = tensorProduct(Nlr(l, rl), nlr(l, rl));
                Ajr(j, R) += rho_aL * M + rho_aT * (I - M);
              }
            }
          });
    
        return Ajr;
      }
    
      NodeValuePerCell<const Rdxd>
      _computeAjr(const SolverType& solver_type,
                  const MeshType& mesh,
                  const DiscreteScalarFunction& rhoaL,
                  const DiscreteScalarFunction& rhoaT) const
      {
        if constexpr (Dimension == 1) {
          return _computeGlaceAjr(mesh, rhoaL, rhoaT);
        } else {
          switch (solver_type) {
          case SolverType::Glace: {
            return _computeGlaceAjr(mesh, rhoaL, rhoaT);
          }
          case SolverType::Eucclhyd: {
            return _computeEucclhydAjr(mesh, rhoaL, rhoaT);
          }
          default: {
            throw UnexpectedError("invalid solver type");
          }
          }
        }
      }
    
      NodeValue<Rdxd>
      _computeAr(const MeshType& mesh, const NodeValuePerCell<const Rdxd>& Ajr) const
      {
        const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
        const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
    
        NodeValue<Rdxd> Ar{mesh.connectivity()};
    
        parallel_for(
          mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
            Rdxd sum                                   = zero;
            const auto& node_to_cell                   = node_to_cell_matrix[r];
            const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
    
            for (size_t j = 0; j < node_to_cell.size(); ++j) {
              const CellId J       = node_to_cell[j];
              const unsigned int R = node_local_number_in_its_cells[j];
              sum += Ajr(J, R);
            }
            Ar[r] = sum;
          });
    
        return Ar;
      }
    
      NodeValue<Rd>
      _computeBr(const Mesh<Connectivity<Dimension>>& mesh,
                 const NodeValuePerCell<const Rdxd>& Ajr,
                 const DiscreteVectorFunction& u,
                 const DiscreteTensorFunction& sigma) const
      {
        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
    
        const NodeValuePerCell<const Rd>& Cjr = mesh_data.Cjr();
    
        const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
        const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
    
        NodeValue<Rd> b{mesh.connectivity()};
    
        parallel_for(
          mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
            const auto& node_to_cell                   = node_to_cell_matrix[r];
            const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(r);
    
            Rd br = zero;
            for (size_t j = 0; j < node_to_cell.size(); ++j) {
              const CellId J       = node_to_cell[j];
              const unsigned int R = node_local_number_in_its_cells[j];
              br += Ajr(J, R) * u[J] - sigma[J] * Cjr(J, R);
            }
    
            b[r] = br;
          });
    
        return b;
      }
    
      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) {
          bool is_valid_boundary_condition = true;
    
          switch (bc_descriptor->type()) {
          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;
          }
          case IBoundaryConditionDescriptor::Type::dirichlet: {
            const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
              dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
            if (dirichlet_bc_descriptor.name() == "velocity") {
              MeshNodeBoundary<Dimension> mesh_node_boundary =
                getMeshNodeBoundary(mesh, dirichlet_bc_descriptor.boundaryDescriptor());
    
              Array<const Rd> value_list =
                InterpolateItemValue<Rd(Rd)>::template interpolate<ItemType::node>(dirichlet_bc_descriptor.rhsSymbolId(),
                                                                                   mesh.xr(),
                                                                                   mesh_node_boundary.nodeList());
    
              bc_list.emplace_back(VelocityBoundaryCondition{mesh_node_boundary, value_list});
            } else if (dirichlet_bc_descriptor.name() == "pressure") {
              const FunctionSymbolId pressure_id = dirichlet_bc_descriptor.rhsSymbolId();
    
              if constexpr (Dimension == 1) {
                MeshNodeBoundary<Dimension> mesh_node_boundary =
                  getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor());
    
                Array<const double> node_values =
                  InterpolateItemValue<double(Rd)>::template interpolate<ItemType::node>(pressure_id, mesh.xr(),
                                                                                         mesh_node_boundary.nodeList());
    
                bc_list.emplace_back(PressureBoundaryCondition{mesh_node_boundary, node_values});
              } else {
                MeshFaceBoundary<Dimension> mesh_face_boundary =
                  getMeshFaceBoundary(mesh, bc_descriptor->boundaryDescriptor());
    
                MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
                Array<const double> face_values =
                  InterpolateItemValue<double(Rd)>::template interpolate<ItemType::face>(pressure_id, mesh_data.xl(),
                                                                                         mesh_face_boundary.faceList());
                bc_list.emplace_back(PressureBoundaryCondition{mesh_face_boundary, face_values});
              }
    
            } else {
              is_valid_boundary_condition = false;
            }
            break;
          }
          default: {
            is_valid_boundary_condition = false;
          }
          }
          if (not is_valid_boundary_condition) {
            std::ostringstream error_msg;
            error_msg << *bc_descriptor << " is an invalid boundary condition for hyperelastic solver";
            throw NormalError(error_msg.str());
          }
        }
    
        return bc_list;
      }
    
      void _applyPressureBC(const BoundaryConditionList& bc_list, const MeshType& mesh, NodeValue<Rd>& br) const;
      void _applySymmetryBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
      void _applyVelocityBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
      void
      _applyBoundaryConditions(const BoundaryConditionList& bc_list,
                               const MeshType& mesh,
                               NodeValue<Rdxd>& Ar,
                               NodeValue<Rd>& br) const
      {
        this->_applyPressureBC(bc_list, mesh, br);
        this->_applySymmetryBC(bc_list, Ar, br);
        this->_applyVelocityBC(bc_list, Ar, br);
      }
    
      NodeValue<const Rd>
      _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
      {
        NodeValue<Rd> u{mesh.connectivity()};
        parallel_for(
          mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { u[r] = inverse(Ar[r]) * br[r]; });
    
        return u;
      }
    
      NodeValuePerCell<Rd>
      _computeFjr(const MeshType& mesh,
                  const NodeValuePerCell<const Rdxd>& Ajr,
                  const NodeValue<const Rd>& ur,
                  const DiscreteVectorFunction& u,
                  const DiscreteTensorFunction& sigma) const
      {
        MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(mesh);
    
        const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    
        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
    
        NodeValuePerCell<Rd> F{mesh.connectivity()};
        parallel_for(
          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
            const auto& cell_nodes = cell_to_node_matrix[j];
    
            for (size_t r = 0; r < cell_nodes.size(); ++r) {
              F(j, r) = -Ajr(j, r) * (u[j] - ur[cell_nodes[r]]) + sigma[j] * Cjr(j, r);
            }
          });
    
        return F;
      }
    
     public:
      std::tuple<const std::shared_ptr<const ItemValueVariant>, const std::shared_ptr<const SubItemValuePerItemVariant>>
      compute_fluxes(const SolverType& solver_type,
                     const std::shared_ptr<const IDiscreteFunction>& i_rho,
                     const std::shared_ptr<const IDiscreteFunction>& i_aL,
                     const std::shared_ptr<const IDiscreteFunction>& i_aT,
                     const std::shared_ptr<const IDiscreteFunction>& i_u,
                     const std::shared_ptr<const IDiscreteFunction>& i_sigma,
                     const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
      {
        std::shared_ptr i_mesh = getCommonMesh({i_rho, i_aL, i_aT, i_u, i_sigma});
        if (not i_mesh) {
          throw NormalError("discrete functions are not defined on the same mesh");
        }
    
        if (not checkDiscretizationType({i_rho, i_aL, i_u, i_sigma}, DiscreteFunctionType::P0)) {
          throw NormalError("hyperelastic solver expects P0 functions");
        }
    
        const MeshType& mesh                = dynamic_cast<const MeshType&>(*i_mesh);
        const DiscreteScalarFunction& rho   = dynamic_cast<const DiscreteScalarFunction&>(*i_rho);
        const DiscreteScalarFunction& aL    = dynamic_cast<const DiscreteScalarFunction&>(*i_aL);
        const DiscreteScalarFunction& aT    = dynamic_cast<const DiscreteScalarFunction&>(*i_aT);
        const DiscreteVectorFunction& u     = dynamic_cast<const DiscreteVectorFunction&>(*i_u);
        const DiscreteTensorFunction& sigma = dynamic_cast<const DiscreteTensorFunction&>(*i_sigma);
    
        NodeValuePerCell<const Rdxd> Ajr = this->_computeAjr(solver_type, mesh, rho * aL, rho * aT);
    
        NodeValue<Rdxd> Ar = this->_computeAr(mesh, Ajr);
        NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, sigma);
    
        const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
        this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
    
        synchronize(Ar);
        synchronize(br);
    
        NodeValue<const Rd> ur         = this->_computeUr(mesh, Ar, br);
        NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
    
        return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
                               std::make_shared<const SubItemValuePerItemVariant>(Fjr));
      }
    
      std::tuple<std::shared_ptr<const IMesh>,
                 std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
                 std::shared_ptr<const DiscreteFunctionP0<Dimension, Rd>>,
                 std::shared_ptr<const DiscreteFunctionP0<Dimension, double>>,
                 std::shared_ptr<const DiscreteFunctionP0<Dimension, Rdxd>>>
      apply_fluxes(const double& dt,
                   const MeshType& mesh,
                   const DiscreteFunctionP0<Dimension, double>& rho,
                   const DiscreteFunctionP0<Dimension, Rd>& u,
                   const DiscreteFunctionP0<Dimension, double>& E,
                   const DiscreteFunctionP0<Dimension, Rdxd>& CG,
                   const NodeValue<const Rd>& ur,
                   const NodeValuePerCell<const Rd>& Fjr) const
      {
        const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
    
        if ((mesh.shared_connectivity() != ur.connectivity_ptr()) or
            (mesh.shared_connectivity() != Fjr.connectivity_ptr())) {
          throw NormalError("fluxes are not defined on the same connectivity than the mesh");
        }
    
        NodeValue<Rd> new_xr = copy(mesh.xr());
        parallel_for(
          mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
    
        std::shared_ptr<const MeshType> new_mesh = std::make_shared<MeshType>(mesh.shared_connectivity(), new_xr);
    
        CellValue<const double> Vj           = MeshDataManager::instance().getMeshData(mesh).Vj();
        const NodeValuePerCell<const Rd> Cjr = MeshDataManager::instance().getMeshData(mesh).Cjr();
    
        CellValue<double> new_rho = copy(rho.cellValues());
        CellValue<Rd> new_u       = copy(u.cellValues());
        CellValue<double> new_E   = copy(E.cellValues());
        CellValue<Rdxd> new_CG    = copy(CG.cellValues());
    
        parallel_for(
          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) {
            const auto& cell_nodes = cell_to_node_matrix[j];
    
            Rd momentum_fluxes       = zero;
            double energy_fluxes     = 0;
            Rdxd cauchy_green_fluxes = zero;
            for (size_t R = 0; R < cell_nodes.size(); ++R) {
              const NodeId r   = cell_nodes[R];
              const Rdxd gradv = tensorProduct(ur[r], Cjr(j, R));
              momentum_fluxes += Fjr(j, R);
              energy_fluxes += dot(Fjr(j, R), ur[r]);
              cauchy_green_fluxes += gradv * CG[j] + CG[j] * transpose(gradv);
            }
            const double dt_over_Mj = dt / (rho[j] * Vj[j]);
            new_u[j] += dt_over_Mj * momentum_fluxes;
            new_E[j] += dt_over_Mj * energy_fluxes;
            new_CG[j] += dt_over_Mj * cauchy_green_fluxes;
            new_CG[j] += transpose(new_CG[j]);
            new_CG[j] *= 0.5;
          });
    
        CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
    
        parallel_for(
          mesh.numberOfCells(), PUGS_LAMBDA(CellId j) { new_rho[j] *= Vj[j] / new_Vj[j]; });
    
        return {new_mesh, std::make_shared<DiscreteScalarFunction>(new_mesh, new_rho),
                std::make_shared<DiscreteVectorFunction>(new_mesh, new_u),
                std::make_shared<DiscreteScalarFunction>(new_mesh, new_E),
                std::make_shared<DiscreteTensorFunction>(new_mesh, new_CG)};
      }
    
      std::tuple<std::shared_ptr<const IMesh>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>>
      apply_fluxes(const double& dt,
                   const std::shared_ptr<const IDiscreteFunction>& rho,
                   const std::shared_ptr<const IDiscreteFunction>& u,
                   const std::shared_ptr<const IDiscreteFunction>& E,
                   const std::shared_ptr<const IDiscreteFunction>& CG,
                   const std::shared_ptr<const ItemValueVariant>& ur,
                   const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr) const
      {
        std::shared_ptr i_mesh = getCommonMesh({rho, u, E});
        if (not i_mesh) {
          throw NormalError("discrete functions are not defined on the same mesh");
        }
    
        if (not checkDiscretizationType({rho, u, E}, DiscreteFunctionType::P0)) {
          throw NormalError("hyperelastic solver expects P0 functions");
        }
    
        return this->apply_fluxes(dt,                                                  //
                                  dynamic_cast<const MeshType&>(*i_mesh),              //
                                  dynamic_cast<const DiscreteScalarFunction&>(*rho),   //
                                  dynamic_cast<const DiscreteVectorFunction&>(*u),     //
                                  dynamic_cast<const DiscreteScalarFunction&>(*E),     //
                                  dynamic_cast<const DiscreteTensorFunction&>(*CG),    //
                                  ur->get<NodeValue<const Rd>>(),                      //
                                  Fjr->get<NodeValuePerCell<const Rd>>());
      }
    
      std::tuple<std::shared_ptr<const IMesh>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>,
                 std::shared_ptr<const IDiscreteFunction>>
      apply(const SolverType& solver_type,
            const double& dt,
            const std::shared_ptr<const IDiscreteFunction>& rho,
            const std::shared_ptr<const IDiscreteFunction>& u,
            const std::shared_ptr<const IDiscreteFunction>& E,
            const std::shared_ptr<const IDiscreteFunction>& CG,
            const std::shared_ptr<const IDiscreteFunction>& aL,
            const std::shared_ptr<const IDiscreteFunction>& aT,
            const std::shared_ptr<const IDiscreteFunction>& sigma,
            const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const
      {
        auto [ur, Fjr] = compute_fluxes(solver_type, rho, aL, aT, u, sigma, bc_descriptor_list);
        return apply_fluxes(dt, rho, u, E, CG, ur, Fjr);
      }
    
      HyperelasticSolver()                     = default;
      HyperelasticSolver(HyperelasticSolver&&) = default;
      ~HyperelasticSolver()                    = default;
    };
    
    template <size_t Dimension>
    void
    HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyPressureBC(const BoundaryConditionList& bc_list,
                                                                               const MeshType& mesh,
                                                                               NodeValue<Rd>& br) const
    {
      for (const auto& boundary_condition : bc_list) {
        std::visit(
          [&](auto&& bc) {
            using T = std::decay_t<decltype(bc)>;
            if constexpr (std::is_same_v<PressureBoundaryCondition, T>) {
              MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(mesh);
              if constexpr (Dimension == 1) {
                const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
    
                const auto& node_to_cell_matrix               = mesh.connectivity().nodeToCellMatrix();
                const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
    
                const auto& node_list  = bc.nodeList();
                const auto& value_list = bc.valueList();
                parallel_for(
                  node_list.size(), PUGS_LAMBDA(size_t i_node) {
                    const NodeId node_id       = node_list[i_node];
                    const auto& node_cell_list = node_to_cell_matrix[node_id];
                    Assert(node_cell_list.size() == 1);
    
                    CellId node_cell_id              = node_cell_list[0];
                    size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
    
                    br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
                  });
              } else {
                const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
    
                const auto& face_to_cell_matrix               = mesh.connectivity().faceToCellMatrix();
                const auto& face_to_node_matrix               = mesh.connectivity().faceToNodeMatrix();
                const auto& face_local_numbers_in_their_cells = mesh.connectivity().faceLocalNumbersInTheirCells();
                const auto& face_cell_is_reversed             = mesh.connectivity().cellFaceIsReversed();
    
                const auto& face_list  = bc.faceList();
                const auto& value_list = bc.valueList();
                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
                  const FaceId face_id       = face_list[i_face];
                  const auto& face_cell_list = face_to_cell_matrix[face_id];
                  Assert(face_cell_list.size() == 1);
    
                  CellId face_cell_id              = face_cell_list[0];
                  size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
    
                  const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
    
                  const auto& face_nodes = face_to_node_matrix[face_id];
    
                  for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
                    NodeId node_id = face_nodes[i_node];
                    br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
                  }
                }
              }
            }
          },
          boundary_condition);
      }
    }
    
    template <size_t Dimension>
    void
    HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applySymmetryBC(const BoundaryConditionList& bc_list,
                                                                               NodeValue<Rdxd>& Ar,
                                                                               NodeValue<Rd>& br) const
    {
      for (const auto& boundary_condition : bc_list) {
        std::visit(
          [&](auto&& bc) {
            using T = std::decay_t<decltype(bc)>;
            if constexpr (std::is_same_v<SymmetryBoundaryCondition, T>) {
              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(
                bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
                  const NodeId r = node_list[r_number];
    
                  Ar[r] = P * Ar[r] * P + nxn;
                  br[r] = P * br[r];
                });
            }
          },
          boundary_condition);
      }
    }
    
    template <size_t Dimension>
    void
    HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyVelocityBC(const BoundaryConditionList& bc_list,
                                                                               NodeValue<Rdxd>& Ar,
                                                                               NodeValue<Rd>& br) const
    {
      for (const auto& boundary_condition : bc_list) {
        std::visit(
          [&](auto&& bc) {
            using T = std::decay_t<decltype(bc)>;
            if constexpr (std::is_same_v<VelocityBoundaryCondition, T>) {
              const auto& node_list  = bc.nodeList();
              const auto& value_list = bc.valueList();
    
              parallel_for(
                node_list.size(), PUGS_LAMBDA(size_t i_node) {
                  NodeId node_id    = node_list[i_node];
                  const auto& value = value_list[i_node];
    
                  Ar[node_id] = identity;
                  br[node_id] = value;
                });
            } else if constexpr (std::is_same_v<FixedBoundaryCondition, T>) {
              const auto& node_list = bc.nodeList();
              parallel_for(
                node_list.size(), PUGS_LAMBDA(size_t i_node) {
                  NodeId node_id = node_list[i_node];
    
                  Ar[node_id] = identity;
                  br[node_id] = zero;
                });
            }
          },
          boundary_condition);
      }
    }
    
    template <size_t Dimension>
    class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::FixedBoundaryCondition
    {
     private:
      const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
    
     public:
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_node_boundary.nodeList();
      }
    
      FixedBoundaryCondition(const MeshNodeBoundary<Dimension> mesh_node_boundary)
        : m_mesh_node_boundary{mesh_node_boundary}
      {}
    
      ~FixedBoundaryCondition() = default;
    };
    
    template <size_t Dimension>
    class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::PressureBoundaryCondition
    {
     private:
      const MeshFaceBoundary<Dimension> m_mesh_face_boundary;
      const Array<const double> m_value_list;
    
     public:
      const Array<const FaceId>&
      faceList() const
      {
        return m_mesh_face_boundary.faceList();
      }
    
      const Array<const double>&
      valueList() const
      {
        return m_value_list;
      }
    
      PressureBoundaryCondition(const MeshFaceBoundary<Dimension>& mesh_face_boundary,
                                const Array<const double>& value_list)
        : m_mesh_face_boundary{mesh_face_boundary}, m_value_list{value_list}
      {}
    
      ~PressureBoundaryCondition() = default;
    };
    
    template <>
    class HyperelasticSolverHandler::HyperelasticSolver<1>::PressureBoundaryCondition
    {
     private:
      const MeshNodeBoundary<1> m_mesh_node_boundary;
      const Array<const double> m_value_list;
    
     public:
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_node_boundary.nodeList();
      }
    
      const Array<const double>&
      valueList() const
      {
        return m_value_list;
      }
    
      PressureBoundaryCondition(const MeshNodeBoundary<1>& mesh_node_boundary, const Array<const double>& value_list)
        : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
      {}
    
      ~PressureBoundaryCondition() = default;
    };
    
    template <size_t Dimension>
    class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::VelocityBoundaryCondition
    {
     private:
      const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
    
      const Array<const TinyVector<Dimension>> m_value_list;
    
     public:
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_node_boundary.nodeList();
      }
    
      const Array<const TinyVector<Dimension>>&
      valueList() const
      {
        return m_value_list;
      }
    
      VelocityBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary,
                                const Array<const TinyVector<Dimension>>& value_list)
        : m_mesh_node_boundary{mesh_node_boundary}, m_value_list{value_list}
      {}
    
      ~VelocityBoundaryCondition() = default;
    };
    
    template <size_t Dimension>
    class HyperelasticSolverHandler::HyperelasticSolver<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();
      }
    
      size_t
      numberOfNodes() const
      {
        return m_mesh_flat_node_boundary.nodeList().size();
      }
    
      const Array<const NodeId>&
      nodeList() const
      {
        return m_mesh_flat_node_boundary.nodeList();
      }
    
      SymmetryBoundaryCondition(const MeshFlatNodeBoundary<Dimension>& mesh_flat_node_boundary)
        : m_mesh_flat_node_boundary(mesh_flat_node_boundary)
      {
        ;
      }
    
      ~SymmetryBoundaryCondition() = default;
    };
    
    HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
    {
      if (not i_mesh) {
        throw NormalError("discrete functions are not defined on the same mesh");
      }
    
      switch (i_mesh->dimension()) {
      case 1: {
        m_hyperelastic_solver = std::make_unique<HyperelasticSolver<1>>();
        break;
      }
      case 2: {
        m_hyperelastic_solver = std::make_unique<HyperelasticSolver<2>>();
        break;
      }
      case 3: {
        m_hyperelastic_solver = std::make_unique<HyperelasticSolver<3>>();
        break;
      }
      default: {
        throw UnexpectedError("invalid mesh dimension");
      }
      }
    }