Add new coupling BC for hyperelastic solver

src/language/modules/SchemeModule.cpp

@@ -815,6 +815,49 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_dt_hyperelastic_eucclhyd_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Eucclhyd, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
                             std::function(
 
@@ -858,6 +901,49 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("local_dt_hyperelastic_glace_solver",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list1,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list2,
+                                 const double& mu, const double& lambda, const double& dt1)
+                                -> std::tuple<
+                                  std::shared_ptr<const MeshVariant>, std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>, std::shared_ptr<const MeshVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>,
+                                  std::shared_ptr<const DiscreteFunctionVariant>> {
+                                return LocalDtHyperelasticSolverHandler{getCommonMesh(
+                                                                          {rho1, u1, E1, CG1, aL1, aT1, sigma1}),
+                                                                        getCommonMesh(
+                                                                          {rho2, u2, E2, CG2, aL2, aT2, sigma2})}
+                                  .solver()
+                                  .apply(LocalDtHyperelasticSolverHandler::SolverType::Glace, dt1, rho1, rho2, u1,
+                                         u2, E1, E2, CG1, CG2, aL1, aL2, aT1, aT2, sigma1, sigma2, bc_descriptor_list1,
+                                         bc_descriptor_list2, mu, lambda);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("apply_hyperelastic_fluxes",
                             std::function(
 

src/scheme/LocalDtHyperelasticSolver.cpp

@@ -16,6 +16,7 @@
 #include <scheme/DiscreteFunctionVariant.hpp>
 #include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
+#include <scheme/HyperelasticSolver.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
@@ -328,6 +329,12 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
                         NodeValuePerCell<Rd>& Fjr,
                         NodeValuePerCell<Rd>& CR_Fjr,
                         const std::vector<NodeId>& map2) const;
+  void _applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+                         NodeValue<Rdxd>& Ar2,
+                         NodeValue<Rd>& ur1,
+                         NodeValue<Rd>& ur2,
+                         const std::vector<NodeId>& map1,
+                         const std::vector<NodeId>& map2) const;
   void
   _applyBoundaryConditions(const BoundaryConditionList& bc_list,
                            const MeshType& mesh,
@@ -608,6 +615,91 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
                            std::make_shared<const SubItemValuePerItemVariant>(Fjr));
   }
 
+    std::tuple<const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>,
+             const std::shared_ptr<const ItemValueVariant>,
+             const std::shared_ptr<const SubItemValuePerItemVariant>,
+             NodeValue<Rd>,
+             NodeValuePerCell<Rd>>
+  compute_fluxes2(const SolverType& solver_type,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u1_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma1_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& rho2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aL2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& aT2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& u2_v,
+                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma2_v,
+                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+                 const std::vector<NodeId>& map1,
+                 const std::vector<NodeId>& map2) const
+  {
+    std::shared_ptr i_mesh1 = getCommonMesh({rho1_v, aL1_v, aT1_v, u1_v, sigma1_v});
+    std::shared_ptr i_mesh2 = getCommonMesh({rho2_v, aL2_v, aT2_v, u2_v, sigma2_v});
+    if (not i_mesh1) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho1_v, aL1_v, u1_v, sigma1_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("hyperelastic solver expects P0 functions");
+    }
+    if (not i_mesh2) {
+      throw NormalError("discrete functions are not defined on the same mesh");
+    }
+
+    if (not checkDiscretizationType({rho2_v, aL2_v, u2_v, sigma2_v}, DiscreteFunctionType::P0)) {
+      throw NormalError("acoustic solver expects P0 functions");
+    }
+
+    const MeshType& mesh1                = *i_mesh1->get<MeshType>();
+    const DiscreteScalarFunction& rho1   = rho1_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u1     = u1_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL1    = aL1_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT1    = aT1_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma1 = sigma1_v->get<DiscreteTensorFunction>();
+
+    const MeshType& mesh2                = *i_mesh2->get<MeshType>();
+    const DiscreteScalarFunction& rho2   = rho2_v->get<DiscreteScalarFunction>();
+    const DiscreteVectorFunction& u2     = u2_v->get<DiscreteVectorFunction>();
+    const DiscreteScalarFunction& aL2    = aL2_v->get<DiscreteScalarFunction>();
+    const DiscreteScalarFunction& aT2    = aT2_v->get<DiscreteScalarFunction>();
+    const DiscreteTensorFunction& sigma2 = sigma2_v->get<DiscreteTensorFunction>();
+
+    NodeValuePerCell<const Rdxd> Ajr1 = this->_computeAjr(solver_type, mesh1, rho1 * aL1, rho1 * aT1);
+    NodeValuePerCell<const Rdxd> Ajr2 = this->_computeAjr(solver_type, mesh2, rho2 * aL2, rho2 * aT2);
+
+    NodeValue<Rdxd> Ar1 = this->_computeAr(mesh1, Ajr1);
+    NodeValue<Rd> br1   = this->_computeBr(mesh1, Ajr1, u1, sigma1);
+    NodeValue<Rdxd> Ar2 = this->_computeAr(mesh2, Ajr2);
+    NodeValue<Rd> br2   = this->_computeBr(mesh2, Ajr2, u2, sigma2);
+
+    const BoundaryConditionList bc_list1 = this->_getBCList(mesh1, bc_descriptor_list1);
+    const BoundaryConditionList bc_list2 = this->_getBCList(mesh2, bc_descriptor_list2);
+    this->_applyBoundaryConditions(bc_list1, mesh1, Ar1, br1);
+    this->_applyBoundaryConditions(bc_list2, mesh2, Ar2, br2);
+
+    synchronize(Ar1);
+    synchronize(br1);
+    synchronize(Ar2);
+    synchronize(br2);
+
+    NodeValue<Rd> ur1           = this->_computeUr(mesh1, Ar1, br1);
+    NodeValue<Rd> ur2           = this->_computeUr(mesh2, Ar2, br2);
+    this->_applyCouplingBC2(Ar1,Ar2,ur1,ur2,map1,map2);
+    NodeValuePerCell<Rd> Fjr1   = this->_computeFjr(mesh1, Ajr1, ur1, u1, sigma1);
+    NodeValuePerCell<Rd> Fjr2   = this->_computeFjr(mesh2, Ajr2, ur2, u2, sigma2);
+
+    return std::make_tuple(std::make_shared<const ItemValueVariant>(ur1),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr1),
+                           std::make_shared<const ItemValueVariant>(ur2),
+                           std::make_shared<const SubItemValuePerItemVariant>(Fjr2),
+                           ur2,
+                           Fjr2);
+  }
+
   std::tuple<std::shared_ptr<const MeshVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
              std::shared_ptr<const DiscreteFunctionVariant>,
@@ -864,6 +956,110 @@ class LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver final
     return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
   }
 
+  std::tuple<std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const MeshVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>,
+             std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const SolverType& solver_type,
+        const double& dt1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+        const double& mu,
+        const double& lambda) const
+  {
+    std::shared_ptr new_m2 = getCommonMesh({rho2, aL2, aT2, u2, sigma2});
+    std::shared_ptr m1     = getCommonMesh({rho1, aL1, aT1, u1, sigma1});
+
+    std::shared_ptr<const DiscreteFunctionVariant> new_rho2 = rho2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_u2   = u2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_E2   = E2;
+    std::shared_ptr<const DiscreteFunctionVariant> new_CG2  = CG2;
+
+    DiscreteScalarFunction aL_d = aL2->get<DiscreteScalarFunction>();
+    DiscreteScalarFunction aT_d = aT2->get<DiscreteScalarFunction>();
+    const std::shared_ptr<const DiscreteFunctionVariant>& c = std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+
+    double dt2 = 0.4 * hyperelastic_dt(c);
+
+    auto [map1, map2] = _computeMapping(m1, new_m2, bc_descriptor_list1, bc_descriptor_list2);
+
+    auto [ur1, Fjr1, ur2, Fjr2, CR_ur, CR_Fjr] =
+      compute_fluxes2(solver_type, rho1, aL1, aT1, u1, sigma1, bc_descriptor_list1, rho2, aL2, aT2, u2, sigma2,
+                     bc_descriptor_list2, map1, map2);
+    const auto [new_m1, new_rho1, new_u1, new_E1, new_CG1] = apply_fluxes(dt1, rho1, u1, E1, CG1, ur1, Fjr1);
+    std::tie(new_m2,new_rho2,new_u2,new_E2,new_CG2) = apply_fluxes(dt2, rho2, u2, E2, CG2, ur2, Fjr2);
+
+    const double gamma            = 1.4;
+    const TinyMatrix<Dimension> I = identity;
+    double sum_dt                 = dt2;
+
+    size_t fluid = 0;
+    if ((mu == 0) and (lambda == 0)) {
+      fluid = 1;
+    }
+
+    while(sum_dt < dt1) {
+
+      const DiscreteScalarFunction& rho_d = new_rho2->get<DiscreteScalarFunction>();
+      const DiscreteVectorFunction& u_d   = new_u2->get<DiscreteVectorFunction>();
+      const DiscreteScalarFunction& E_d   = new_E2->get<DiscreteScalarFunction>();
+      const DiscreteScalarFunction& eps   = E_d - 0.5 * dot(u_d, u_d);
+      const DiscreteTensorFunction& CG_d  = new_CG2->get<DiscreteTensorFunction>();
+      const DiscreteScalarFunction& p = fluid * (gamma - 1) * rho_d * eps;
+      const DiscreteTensorFunction& sigma_d =
+        mu / sqrt(det(CG_d)) * (CG_d - I) + lambda / sqrt(det(CG_d)) * log(sqrt(det(CG_d))) * I - p * I;
+      aL_d = sqrt((lambda + 2 * mu) / rho_d + gamma * p / rho_d);
+      aT_d = sqrt(mu / rho_d);
+
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_sigma2 =
+        std::make_shared<const DiscreteFunctionVariant>(sigma_d);
+        const std::shared_ptr<const DiscreteFunctionVariant>& new_c =
+        std::make_shared<const DiscreteFunctionVariant>(aL_d + aT_d);
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_aL2 =
+        std::make_shared<const DiscreteFunctionVariant>(aL_d);
+      const std::shared_ptr<const DiscreteFunctionVariant>& new_aT2 =
+        std::make_shared<const DiscreteFunctionVariant>(aT_d);
+
+      dt2 = 0.4 * hyperelastic_dt(new_c);
+      if(sum_dt + dt2 > dt1){
+	    dt2 = dt1 - sum_dt;
+      }
+
+      auto [sub_ur2, sub_Fjr2] = compute_fluxes(solver_type, new_rho2, new_aL2, new_aT2, new_u2, new_sigma2,
+                                        bc_descriptor_list2, CR_ur, CR_Fjr, map2);
+
+      std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+        apply_fluxes(dt2, new_rho2, new_u2, new_E2, new_CG2, sub_ur2, sub_Fjr2);
+
+      sum_dt += dt2;
+    }
+
+    std::tie(new_m2, new_rho2, new_u2, new_E2, new_CG2) =
+      mesh_correction(new_m1, new_m2, new_rho2, new_u2, new_E2, new_CG2, map1, map2);
+
+    return {new_m1, new_rho1, new_u1, new_E1, new_CG1, new_m2, new_rho2, new_u2, new_E2, new_CG2};
+  }
+
   LocalDtHyperelasticSolver()                            = default;
   LocalDtHyperelasticSolver(LocalDtHyperelasticSolver&&) = default;
   ~LocalDtHyperelasticSolver()                           = default;
@@ -1095,6 +1291,30 @@ LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCou
   }
 }
 
+template <MeshConcept MeshType>
+void 
+LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC2(NodeValue<Rdxd>& Ar1,
+										      NodeValue<Rdxd>& Ar2,
+										      NodeValue<Rd>& ur1,
+										      NodeValue<Rd>& ur2,
+										      const std::vector<NodeId>& map1,
+										      const std::vector<NodeId>& map2) const
+{
+  size_t n = map1.size();
+  Rd lambda;
+
+  for(size_t i = 0; i<n; i++){
+
+    NodeId node_id1 = map1[i];
+    NodeId node_id2 = map2[i];
+    
+    lambda = inverse(inverse(Ar2[node_id2]) + inverse(Ar1[node_id1]))*(ur2[node_id2]-ur1[node_id1]); 
+
+    ur1[node_id1] = ur1[node_id1] + inverse(Ar1[node_id1])*lambda;
+    ur2[node_id2] = ur2[node_id2] - inverse(Ar2[node_id2])*lambda; 
+  }
+}
+
 template <MeshConcept MeshType>
 void
 LocalDtHyperelasticSolverHandler::LocalDtHyperelasticSolver<MeshType>::_applyCouplingBC(

src/scheme/LocalDtHyperelasticSolver.hpp

@@ -82,6 +82,37 @@ class LocalDtHyperelasticSolverHandler
           const double& mu,
           const double& lambda) const = 0;
 
+    virtual std::tuple<std::shared_ptr<const MeshVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const MeshVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>,
+                       std::shared_ptr<const DiscreteFunctionVariant>>
+  apply(const SolverType& solver_type,
+        const double& dt1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& rho2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& u2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& E2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& CG2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aL2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& aT2,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma1,
+        const std::shared_ptr<const DiscreteFunctionVariant>& sigma2,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list1,
+        const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list2,
+        const double& mu,
+        const double& lambda) const = 0;
+
     ILocalDtHyperelasticSolver()                                        = default;
     ILocalDtHyperelasticSolver(ILocalDtHyperelasticSolver&&)            = default;
     ILocalDtHyperelasticSolver& operator=(ILocalDtHyperelasticSolver&&) = default;