diff --git a/src/scheme/HyperelasticSolver.cpp b/src/scheme/HyperelasticSolver.cpp
index c18ca1e54be07e11b8b94f1a03efb6da6b58af37..dab77ca4cb74d692f3e63bd03cf85b18729271a6 100644
--- a/src/scheme/HyperelasticSolver.cpp
+++ b/src/scheme/HyperelasticSolver.cpp
@@ -421,6 +421,10 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
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 _applyCouplingBC2(const MeshType& mesh,
+ const BoundaryConditionList& bc_list,
+ NodeValue<Rd>& ur,
+ NodeValuePerCell<Rd>& Fjr) const;
void
_applyBoundaryConditions(const BoundaryConditionList& bc_list,
const MeshType& mesh,
@@ -436,7 +440,7 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
this->_applyVelocityBC(bc_list, Ar, br);
}
- NodeValue<const Rd>
+ NodeValue<Rd>
_computeUr(const MeshType& mesh, const NodeValue<Rdxd>& Ar, const NodeValue<Rd>& br) const
{
NodeValue<Rd> u{mesh.connectivity()};
@@ -510,8 +514,10 @@ class HyperelasticSolverHandler::HyperelasticSolver final : public HyperelasticS
synchronize(Ar);
synchronize(br);
- NodeValue<const Rd> ur = this->_computeUr(mesh, Ar, br);
- NodeValuePerCell<const Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+ NodeValue<Rd> ur = this->_computeUr(mesh, Ar, br);
+ NodeValuePerCell<Rd> Fjr = this->_computeFjr(mesh, Ajr, ur, u, sigma);
+
+ this->_applyCouplingBC2(mesh, bc_list, ur, Fjr);
return std::make_tuple(std::make_shared<const ItemValueVariant>(ur),
std::make_shared<const SubItemValuePerItemVariant>(Fjr));
@@ -946,6 +952,141 @@ HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC(const
}
}
+template <size_t Dimension>
+void
+HyperelasticSolverHandler::HyperelasticSolver<Dimension>::_applyCouplingBC2(const MeshType& mesh,
+ const BoundaryConditionList& bc_list,
+ NodeValue<Rd>& ur,
+ NodeValuePerCell<Rd>& Fjr) 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;
+ auto rank = costo->myGlobalRank();
+ if (rank == 0) {
+ const size_t& n = node_list.size();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ for (size_t i = 0; i < n; i++) {
+ const NodeId node_id = node_list[i];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+ /* std::cout << "\033[01;33m" */
+ /* << "rank_0: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
+ 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];
+ /* Fjr(J, R) = CR_Fjr(J, R); */
+ }
+ /* ur[node_id] = CR_ur[node_id]; */
+ }
+ }
+ if (rank == 1) {
+ /* std::vector<double> recv; */
+ /* tag = 310; */
+ /* costo->recvData(recv, dest, tag); */
+
+ /* CR_Fjr = ... */
+ /* CR_ur = */
+ const size_t& n = node_list.size();
+
+ const auto& node_to_cell_matrix = mesh.connectivity().nodeToCellMatrix();
+ const auto& node_local_numbers_in_their_cells = mesh.connectivity().nodeLocalNumbersInTheirCells();
+
+ for (size_t i = 0; i < n; i++) {
+ const NodeId node_id = node_list[i];
+ const auto& node_to_cell = node_to_cell_matrix[node_id];
+ const auto& node_local_number_in_its_cells = node_local_numbers_in_their_cells.itemArray(node_id);
+ /* std::cout << "\033[01;32m" */
+ /* << "rank_1: " << node_id << " " << node_to_cell.size() << "\033[00;00m" << std::endl; */
+
+ 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];
+ /* Fjr(J, R) = CR_Fjr(J, R); */
+ }
+ /* ur[node_id] = CR_ur[node_id]; */
+ }
+ }
+
+ /* const int dest = costo->myGlobalSize() - 1; */
+ /* int tag = 300; */
+ /* 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 = 310; */
+ /* 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
{