From 41124c36ebe1b599d685b7c085cd9eb767172077 Mon Sep 17 00:00:00 2001
From: "t. chantrait" <teddy.chantrait@cea.fr>
Date: Mon, 2 Jun 2025 13:52:07 +0200
Subject: [PATCH] Fix symmetric BC
---
src/language/modules/SchemeModule.cpp | 47 ++++++++++++++-----
src/scheme/SerrailleSolver.cpp | 65 ++++++++++++++++++++++-----
src/scheme/SerrailleSolver.hpp | 4 +-
3 files changed, 94 insertions(+), 22 deletions(-)
diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 62744a6d9..fafb9a4b3 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/src/language/modules/SchemeModule.cpp
@@ -1531,7 +1531,7 @@ SchemeModule::SchemeModule()
));
- this->_addBuiltinFunction("hyperplastic_glace_solver_serraille",
+ this->_addBuiltinFunction("apply_hyperplastic_fluxes_cauchygreen_serraille",
std::function(
[](const std::shared_ptr<const DiscreteFunctionVariant>& rho,
@@ -1540,9 +1540,9 @@ SchemeModule::SchemeModule()
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>& aL,
- const std::shared_ptr<const DiscreteFunctionVariant>& aT,
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>,
@@ -1550,15 +1550,40 @@ SchemeModule::SchemeModule()
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>,
std::shared_ptr<const DiscreteFunctionVariant>> {
- return HyperplasticSolverSerrailleHandler{
- getCommonMesh({rho, u, E, Fe, zeta, yield, aL, aT, sigma})}
- .solver()
- .apply(HyperplasticSolverSerrailleHandler::SolverType::Glace,
- HyperplasticSolverSerrailleHandler::RelaxationType::CauchyGreen, dt, rho, u, E, Fe,
- zeta, yield, aL, aT, sigma, bc_descriptor_list);
- }
+ 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(
+
+ [](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>& aL,
+ const std::shared_ptr<const DiscreteFunctionVariant>& aT,
+ const std::shared_ptr<const DiscreteFunctionVariant>& sigma,
+ 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, zeta, yield, aL, aT, sigma})}
+ .solver()
+ .apply(HyperplasticSolverSerrailleHandler::SolverType::Glace,
+ HyperplasticSolverSerrailleHandler::RelaxationType::CauchyGreen, dt, rho, u, E, Fe,
+ zeta, yield, aL, aT, sigma, bc_descriptor_list);
+ }
+
+ ));
diff --git a/src/scheme/SerrailleSolver.cpp b/src/scheme/SerrailleSolver.cpp
index 717a4f093..27912a1ea 100644
--- a/src/scheme/SerrailleSolver.cpp
+++ b/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>,
@@ -1485,6 +1524,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
diff --git a/src/scheme/SerrailleSolver.hpp b/src/scheme/SerrailleSolver.hpp
index 8b61a4c19..a5fb58acd 100644
--- a/src/scheme/SerrailleSolver.hpp
+++ b/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>,
--
GitLab