Hyperelastic Lagrangian solver (Eucclhyd or Glace-type).

First commit.

src/scheme/HyperelasticSolver.cpp

+#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 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;
+      });
+
+    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");
+  }
+  }
+}

src/scheme/HyperelasticSolver.hpp

+#ifndef HYPERELASTIC_SOLVER_HPP
+#define HYPERELASTIC_SOLVER_HPP
+
+#include <memory>
+#include <tuple>
+#include <vector>
+
+class IDiscreteFunction;
+class IBoundaryConditionDescriptor;
+class IMesh;
+class ItemValueVariant;
+class SubItemValuePerItemVariant;
+
+double hyperelastic_dt(const std::shared_ptr<const IDiscreteFunction>& c);
+
+class HyperelasticSolverHandler
+{
+ public:
+  enum class SolverType
+  {
+    Glace,
+    Eucclhyd
+  };
+
+ private:
+  struct IHyperelasticSolver
+  {
+    virtual 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>& rho,
+      const std::shared_ptr<const IDiscreteFunction>& aL,
+      const std::shared_ptr<const IDiscreteFunction>& aT,
+      const std::shared_ptr<const IDiscreteFunction>& u,
+      const std::shared_ptr<const IDiscreteFunction>& sigma,
+      const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+    virtual 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 = 0;
+
+    virtual 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>& p,
+          const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list) const = 0;
+
+    IHyperelasticSolver()                                 = default;
+    IHyperelasticSolver(IHyperelasticSolver&&)            = default;
+    IHyperelasticSolver& operator=(IHyperelasticSolver&&) = default;
+
+    virtual ~IHyperelasticSolver() = default;
+  };
+
+  template <size_t Dimension>
+  class HyperelasticSolver;
+
+  std::unique_ptr<IHyperelasticSolver> m_hyperelastic_solver;
+
+ public:
+  const IHyperelasticSolver&
+  solver() const
+  {
+    return *m_hyperelastic_solver;
+  }
+
+  HyperelasticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
+};
+
+#endif   // HYPERELASTIC_SOLVER_HPP