report in hyperelasticity

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) {

src/scheme/HyperelasticSolver.hpp

-#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