diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index 3a9725d107d801f58c915a090b0ccc840f4bbde8..c18ca1e54be07e11b8b94f1a03efb6da6b58af37 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -7,6 +7,7 @@
#include <mesh/MeshFlatNodeBoundary.hpp>
#include <mesh/MeshNodeBoundary.hpp>
#include <mesh/SubItemValuePerItemVariant.hpp>
+#include <scheme/CouplingBoundaryConditionDescriptor.hpp>
#include <scheme/DirichletBoundaryConditionDescriptor.hpp>
#include <scheme/DiscreteFunctionP0.hpp>
#include <scheme/DiscreteFunctionUtils.hpp>
@@ -81,6 +82,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
#endif // PUGS_HAS_COSTO
class SymmetryBoundaryCondition;
class VelocityBoundaryCondition;
+ class CouplingBoundaryCondition;
using BoundaryCondition = std::variant<FixedBoundaryCondition,
PressureBoundaryCondition,
@@ -88,6 +90,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
#ifdef PUGS_HAS_COSTO
CplDirichletBoundaryCondition,
#endif // PUGS_HAS_COSTO
+ CouplingBoundaryCondition,
SymmetryBoundaryCondition,
VelocityBoundaryCondition>;
@@ -390,6 +393,12 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
}
break;
}
+ case IBoundaryConditionDescriptor::Type::coupling: {
+ const CouplingBoundaryConditionDescriptor& coupling_bc_descriptor =
+ dynamic_cast<const CouplingBoundaryConditionDescriptor&>(*bc_descriptor);
+ bc_list.emplace_back(CouplingBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
+ break;
+ }
default: {
is_valid_boundary_condition = false;
}
@@ -411,6 +420,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
#endif // PUGS_HAS_COSTO
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 _applyCouplingBC(const BoundaryConditionList& bc_list, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br) const;
void
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
const MeshType& mesh,
@@ -495,6 +505,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
const BoundaryConditionList bc_list = this->_getBCList(mesh, bc_descriptor_list);
this->_applyBoundaryConditions(bc_list, mesh, Ar, br);
+ this->_applyCouplingBC(bc_list, Ar, br);
synchronize(Ar);
synchronize(br);
@@ -854,6 +865,87 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyVelocityBC(const
}
}
+template <size_t Dimension>
+void
+HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(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<CouplingBoundaryCondition, T>) {
+ const auto& node_list = bc.nodeList();
+
+ /* std::cout << "\033[01;31m" */
+ /* << "un applyCoupling" */
+ /* << "Dimension: " << Dimension << "\033[00;00m" << std::endl; */
+ /* std::cout << "\033[01;31m" */
+ /* << "node_list.size()" << node_list.size() << "\033[00;00m" << std::endl; */
+
+#ifdef PUGS_HAS_COSTO
+ auto costo = parallel::Messenger::getInstance().myCoupling;
+ const int dest = costo->myGlobalSize() - 1;
+ int tag = 200;
+ std::vector<int> shape;
+ shape.resize(3);
+ shape[0] = node_list.size();
+ shape[1] = Dimension + Dimension * Dimension;
+ /* shape[1] = Dimension; */
+ shape[2] = 0;
+
+ std::vector<double> data;
+ data.resize(shape[0] * shape[1] + shape[2]);
+ Array<TinyVector<Dimension>> br_extracted(node_list.size());
+ Array<TinyMatrix<Dimension>> Ar_extracted(node_list.size());
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br_extracted[i_node] = br[node_id];
+ Ar_extracted[i_node] = Ar[node_id];
+ }
+ });
+ /*TODO: serializer directement Ar et br pour eviter une copie*/
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ data[i_node * shape[1] + i_dim] = br_extracted[i_node][i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ data[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim] = Ar_extracted[i_node](i_dim, j_dim);
+ }
+ }
+ }
+ costo->sendData(shape, data, dest, tag);
+ std::vector<double> recv;
+
+ tag = 210;
+ costo->recvData(recv, dest, tag);
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ for (unsigned short i_dim = 0; i_dim < Dimension; ++i_dim) {
+ br_extracted[i_node][i_dim] = recv[i_node * shape[1] + i_dim];
+ for (size_t j_dim = 0; j_dim < Dimension; j_dim++) {
+ Ar_extracted[i_node](i_dim, j_dim) = recv[i_node * shape[1] + Dimension + i_dim * Dimension + j_dim];
+ }
+ }
+ }
+
+ parallel_for(
+ node_list.size(), PUGS_LAMBDA(size_t i_node) {
+ NodeId node_id = node_list[i_node];
+ for (size_t i_dim = 0; i_dim < Dimension; i_dim++) {
+ br[node_id] = br_extracted[i_node];
+ Ar[node_id] = Ar_extracted[i_node];
+ }
+ });
+
+#endif PUGS_HAS_COSTO
+ }
+ },
+ boundary_condition);
+ }
+}
+
template <size_t Dimension>
class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::FixedBoundaryCondition
{
@@ -1080,6 +1172,28 @@ class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::SymmetryBoundary
~SymmetryBoundaryCondition() = default;
};
+template <size_t Dimension>
+class HyperelasticSolverHandler::HyperelasticSolver<Dimension>::CouplingBoundaryCondition
+{
+ private:
+ const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+
+ public:
+ const Array<const NodeId>&
+ nodeList() const
+ {
+ return m_mesh_node_boundary.nodeList();
+ }
+
+ CouplingBoundaryCondition(const MeshNodeBoundary<Dimension>& mesh_node_boundary)
+ : m_mesh_node_boundary{mesh_node_boundary}
+ {
+ ;
+ }
+
+ ~CouplingBoundaryCondition() = default;
+};
+
HyperelasticSolverHandler::HyperelasticSolverHandler(const std::shared_ptr<const IMesh>& i_mesh)
{
if (not i_mesh) {
diff --git a/src/scheme/HyperelasticSolver.hpp b/src/scheme/HyperelasticSolver.hpp
index f82cbe361051efc64a3e413d60c830febd62d907..044ccd96f80fff4a46613b32346b19b46af23ed1 100644
--- a/src/scheme/HyperelasticSolver.hpp
+++ b/src/scheme/HyperelasticSolver.hpp
@@ -1,6 +1,4 @@
-#ifndef HYPERELASTIC_SOLVER_HPP
-#define HYPERELASTIC_SOLVER_HPP
-
+#pragma once
#include <memory>
#include <tuple>
#include <vector>
@@ -86,5 +84,3 @@ class HyperelasticSolverHandler
HyperelasticSolverHandler(const std::shared_ptr<const IMesh>& mesh);
};
-
-#endif // HYPERELASTIC_SOLVER_HPP