Fix symmetric BC

src/language/modules/SchemeModule.cpp

@@ -1531,6 +1531,31 @@ SchemeModule::SchemeModule()
                               ));
 
 
+  this->_addBuiltinFunction("apply_hyperplastic_fluxes_cauchygreen_serraille",
+                            std::function(
+
+                              [](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& u,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& E,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& Fe,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& zeta,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& yield,
+                                 const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+                                 const std::shared_ptr<const ItemValueVariant>& ur,
+                                 const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
+                                 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>&
+                                   bc_descriptor_list,
+                                 const double& dt) -> 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>> {
+                                return HyperplasticSolverSerrailleHandler{getCommonMesh({rho, u, E, Fe})}                                  .solver()
+                                  .apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, bc_descriptor_list,
+						HyperplasticSolverSerrailleHandler::RelaxationType::CauchyGreen
+						);
+                              }));
+
   this->_addBuiltinFunction("hyperplastic_glace_solver_serraille",
                             std::function(
 

src/scheme/SerrailleSolver.cpp

@@ -434,6 +434,29 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
     this->_applySymmetryBC(bc_list, Ar, br);
     this->_applyVelocityBC(bc_list, Ar, br);
   }
+  void
+  _projectOnBoundary(const BoundaryConditionList& bc_list, NodeValue<Rd>& xr) 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 Rd& o = bc.origin();
+
+            const Array<const NodeId>& node_list = bc.nodeList();
+            parallel_for(
+              bc.numberOfNodes(), PUGS_LAMBDA(int r_number) {
+                const NodeId r = node_list[r_number];
+                Rd& x          = xr[node_list[r]];
+                x -= dot(x - o, n) * n;
+              });
+          }
+        },
+        boundary_condition);
+    }
+  }
 
   NodeValue<const Rd>
   _computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
@@ -536,7 +559,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                const DiscreteScalarFunction& yield,
                const DiscreteTensorFunction3d& sigma,
                const NodeValue<const Rd>& ur,
-               const NodeValuePerCell<const Rd>& Fjr) const
+               const NodeValuePerCell<const Rd>& Fjr,
+	       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list
+	       ) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -623,7 +648,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                 const DiscreteScalarFunction& yield,
                 const DiscreteTensorFunction3d& sigma,
                 const NodeValue<const Rd>& ur,
-                const NodeValuePerCell<const Rd>& Fjr) const
+                const NodeValuePerCell<const Rd>& Fjr,
+		const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list
+		) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -711,7 +738,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                  const DiscreteScalarFunction& yield,
                  const DiscreteTensorFunction3d& sigma,
                  const NodeValue<const Rd>& ur,
-                 const NodeValuePerCell<const Rd>& Fjr) const
+                 const NodeValuePerCell<const Rd>& Fjr,
+		 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list
+		 ) const
   {
     const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
 
@@ -723,6 +752,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
     NodeValue<Rd> new_xr = copy(mesh.xr());
     parallel_for(mesh.numberOfNodes(), PUGS_LAMBDA(NodeId r) { new_xr[r] += dt * ur[r]; });
 
+    const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
+    this->_projectOnBoundary(bc_list, new_xr);
+
     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();
@@ -778,7 +810,6 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
         // const R3x3 expMT = EigenvalueSolver{}.computeExpForTinyMatrix(transpose(M));
 
         // calcul de S trial
-
         const R3x3 Strial   = (new_mu[j] / J) * (Be_barre - 1./3*trace(Be_barre)*I);
         double limit2       = 2. / 3 * yield[j] * yield[j];
         double normS2 = frobeniusNorm(Strial);
@@ -800,7 +831,6 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
         // const double fenorm1 = frobeniusNorm(new_Be[j]);
         // rationorm[j]         = fenorm1 / fenorm0;
       });
-    std::cout << "sum_chi " << sum(chi) << "\n";
 
     CellValue<const double> new_Vj = MeshDataManager::instance().getMeshData(*new_mesh).Vj();
 
@@ -828,7 +858,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                const std::shared_ptr<const ItemValueVariant>& ur,
                const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
-               const RelaxationType& relaxation_type) const
+	       const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,	       
+               const RelaxationType& relaxation_type
+	       ) const
   {
     std::shared_ptr mesh_v = getCommonMesh({rho, u, E});
     if (not mesh_v) {
@@ -850,7 +882,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                                 yield->get<DiscreteScalarFunction>(),     //
                                 sigma->get<DiscreteTensorFunction3d>(),   //
                                 ur->get<NodeValue<const Rd>>(),           //
-                                Fjr->get<NodeValuePerCell<const Rd>>());
+                                Fjr->get<NodeValuePerCell<const Rd>>(),
+				bc_descriptor_list
+				);
     }
     case RelaxationType::Implicit: {
       return this->apply_fluxes2(dt,                                       //
@@ -863,7 +897,9 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                                  yield->get<DiscreteScalarFunction>(),     //
                                  sigma->get<DiscreteTensorFunction3d>(),   //
                                  ur->get<NodeValue<const Rd>>(),           //
-                                 Fjr->get<NodeValuePerCell<const Rd>>());
+                                 Fjr->get<NodeValuePerCell<const Rd>>(),
+				 bc_descriptor_list
+				 );
     }
     case RelaxationType::CauchyGreen: {
       return this->apply_fluxes_B(dt,                                       //
@@ -876,7 +912,10 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
                                   yield->get<DiscreteScalarFunction>(),     //
                                   sigma->get<DiscreteTensorFunction3d>(),   //
                                   ur->get<NodeValue<const Rd>>(),           //
-                                  Fjr->get<NodeValuePerCell<const Rd>>());
+                                  Fjr->get<NodeValuePerCell<const Rd>>(),
+				  bc_descriptor_list
+				  );
+
     }
     default: {
       throw UnexpectedError("invalid relaxation type");
@@ -1088,7 +1127,7 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver final
         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, Fe, zeta, yield, sigma, ur, Fjr, relaxation_type);
+    return apply_fluxes(dt, rho, u, E, Fe, zeta, yield, sigma, ur, Fjr, bc_descriptor_list, relaxation_type);
   }
 
   std::tuple<std::shared_ptr<const MeshVariant>,
@@ -1486,6 +1525,12 @@ class HyperplasticSolverSerrailleHandler::SerrailleSolver<MeshType>::SymmetryBou
     return m_mesh_flat_node_boundary.outgoingNormal();
   }
   
+  const Rd&
+  origin() const
+  {
+    return m_mesh_flat_node_boundary.origin();
+  }
+
   size_t
   numberOfNodes() const
   {

src/scheme/SerrailleSolver.hpp

@@ -59,7 +59,9 @@ class HyperplasticSolverSerrailleHandler
                  const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
                  const std::shared_ptr<const ItemValueVariant>& ur,
                  const std::shared_ptr<const SubItemValuePerItemVariant>& Fjr,
-                 const RelaxationType& relaxation_type) const = 0;
+		 const std::vector<std::shared_ptr<const IBoundaryConditionDescriptor>>& bc_descriptor_list,
+                 const RelaxationType& relaxation_type
+		 ) const = 0;
 
     virtual std::tuple<std::shared_ptr<const MeshVariant>,
                        std::shared_ptr<const DiscreteFunctionVariant>,