Add a generic external bc descriptor

Use the acoustic solver (in 1d) to check the possibility of exchanging
data through sockets. Here, one provides the pressure on a boundary to
an external solver which in returns provides the boundary velocity.

src/language/modules/SchemeModule.cpp

@@ -26,6 +26,7 @@
 #include <scheme/DiscreteFunctionUtils.hpp>
 #include <scheme/DiscreteFunctionVectorIntegrator.hpp>
 #include <scheme/DiscreteFunctionVectorInterpoler.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/FourierBoundaryConditionDescriptor.hpp>
 #include <scheme/FreeBoundaryConditionDescriptor.hpp>
@@ -36,6 +37,7 @@
 #include <scheme/NeumannBoundaryConditionDescriptor.hpp>
 #include <scheme/NumberedBoundaryDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
 
 #include <memory>
 
@@ -273,6 +275,20 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("external_fsi_velocity",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+                                                                  const std::shared_ptr<const Socket>&)>>(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const std::shared_ptr<const Socket>& socket)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<ExternalBoundaryConditionDescriptor>("external_fsi_velocity",
+                                                                                             boundary, socket);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("glace_solver",
                             std::make_shared<BuiltinFunctionEmbedder<std::tuple<
                               std::shared_ptr<const IMesh>, std::shared_ptr<const IDiscreteFunction>,

src/scheme/AcousticSolver.cpp

@@ -8,11 +8,13 @@
 #include <scheme/DirichletBoundaryConditionDescriptor.hpp>
 #include <scheme/DiscreteFunctionP0.hpp>
 #include <scheme/DiscreteFunctionUtils.hpp>
+#include <scheme/ExternalBoundaryConditionDescriptor.hpp>
 #include <scheme/FixedBoundaryConditionDescriptor.hpp>
 #include <scheme/IBoundaryConditionDescriptor.hpp>
 #include <scheme/IDiscreteFunction.hpp>
 #include <scheme/IDiscreteFunctionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
 
 #include <variant>
 #include <vector>
@@ -75,9 +77,13 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   class PressureBoundaryCondition;
   class SymmetryBoundaryCondition;
   class VelocityBoundaryCondition;
+  class ExternalFSIVelocityBoundaryCondition;
 
-  using BoundaryCondition = std::
-    variant<FixedBoundaryCondition, PressureBoundaryCondition, SymmetryBoundaryCondition, VelocityBoundaryCondition>;
+  using BoundaryCondition = std::variant<FixedBoundaryCondition,
+                                         PressureBoundaryCondition,
+                                         SymmetryBoundaryCondition,
+                                         VelocityBoundaryCondition,
+                                         ExternalFSIVelocityBoundaryCondition>;
 
   using BoundaryConditionList = std::vector<BoundaryCondition>;
 
@@ -254,6 +260,18 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
         bc_list.emplace_back(FixedBoundaryCondition(getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor())));
         break;
       }
+      case IBoundaryConditionDescriptor::Type::external: {
+        if constexpr (Dimension == 1) {
+          const ExternalBoundaryConditionDescriptor& extern_bc_descriptor =
+            dynamic_cast<const ExternalBoundaryConditionDescriptor&>(*bc_descriptor);
+          bc_list.emplace_back(
+            ExternalFSIVelocityBoundaryCondition(mesh, getMeshNodeBoundary(mesh, bc_descriptor->boundaryDescriptor()),
+                                                 extern_bc_descriptor.socket()));
+        } else {
+          throw NotImplementedError("external FSI velocity not implemented for dimension > 1");
+        }
+        break;
+      }
       case IBoundaryConditionDescriptor::Type::dirichlet: {
         const DirichletBoundaryConditionDescriptor& dirichlet_bc_descriptor =
           dynamic_cast<const DirichletBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -320,6 +338,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
   void _applyPressureBC(const MeshType& mesh, NodeValue<Rd>& br);
   void _applySymmetryBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
   void _applyVelocityBC(NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
+  void _applyExternalVelocityBC(const DiscreteScalarFunction& p, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br);
 
   void
   _applyBoundaryConditions(const MeshType& mesh, NodeValue<Rdxd>& Ar, NodeValue<Rd>& br)
@@ -380,6 +399,7 @@ class AcousticSolverHandler::AcousticSolver final : public AcousticSolverHandler
     NodeValue<Rd> br   = this->_computeBr(mesh, Ajr, u, p);
 
     this->_applyBoundaryConditions(mesh, Ar, br);
+    this->_applyExternalVelocityBC(p, Ar, br);
 
     synchronize(Ar);
     synchronize(br);
@@ -610,6 +630,34 @@ AcousticSolverHandler::AcousticSolver<Dimension>::_applyVelocityBC(NodeValue<Rdx
   }
 }
 
+template <size_t Dimension>
+void
+AcousticSolverHandler::AcousticSolver<Dimension>::_applyExternalVelocityBC(const DiscreteScalarFunction& p,
+                                                                           NodeValue<Rdxd>& Ar,
+                                                                           NodeValue<Rd>& br)
+{
+  for (const auto& boundary_condition : m_boundary_condition_list) {
+    std::visit(
+      [&](auto&& bc) {
+        using T = std::decay_t<decltype(bc)>;
+        if constexpr (std::is_same_v<ExternalFSIVelocityBoundaryCondition, T>) {
+          const auto& node_list  = bc.nodeList();
+          const auto& value_list = bc.valueList(p);
+
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
+              NodeId node_id    = node_list[i_node];
+              const auto& value = value_list[i_node];
+
+              Ar[node_id] = identity;
+              br[node_id] = value;
+            });
+        }
+      },
+      boundary_condition);
+  }
+}
+
 template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver<Dimension>::FixedBoundaryCondition
 {
@@ -685,6 +733,55 @@ class AcousticSolverHandler::AcousticSolver<1>::PressureBoundaryCondition
   ~PressureBoundaryCondition() = default;
 };
 
+template <size_t Dimension>
+class AcousticSolverHandler::AcousticSolver<Dimension>::ExternalFSIVelocityBoundaryCondition
+{
+ private:
+  const ItemToItemMatrix<ItemType::node, ItemType::cell> m_node_to_cell_matrix;
+  const MeshNodeBoundary<Dimension> m_mesh_node_boundary;
+  const Array<TinyVector<Dimension>> m_value_list;
+  const std::shared_ptr<const Socket> m_socket;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_mesh_node_boundary.nodeList();
+  }
+
+  Array<const TinyVector<Dimension>>
+  valueList(const DiscreteScalarFunction& p) const
+  {
+    if (parallel::size() > 1) {
+      throw NotImplementedError("Parallelism is not supported");
+    }
+    if (m_value_list.size() != 1) {
+      throw NotImplementedError("Non connected boundaries are not supported");
+    }
+
+    const CellId cell_id = m_node_to_cell_matrix[m_mesh_node_boundary.nodeList()[0]][0];
+
+    write(*m_socket, p[cell_id]);
+    read(*m_socket, m_value_list[0]);
+    return m_value_list;
+  }
+
+  ExternalFSIVelocityBoundaryCondition(const Mesh<Connectivity<Dimension>>& mesh,
+                                       const MeshNodeBoundary<Dimension>& mesh_node_boundary,
+                                       const std::shared_ptr<const Socket>& socket)
+    : m_node_to_cell_matrix{mesh.connectivity().nodeToCellMatrix()},
+      m_mesh_node_boundary{mesh_node_boundary},
+      m_value_list(mesh_node_boundary.nodeList().size()),
+      m_socket{socket}
+  {
+    if constexpr (Dimension > 1) {
+      throw NotImplementedError("external FSI velocity bc in dimension > 1");
+    }
+  }
+
+  ~ExternalFSIVelocityBoundaryCondition() = default;
+};
+
 template <size_t Dimension>
 class AcousticSolverHandler::AcousticSolver<Dimension>::VelocityBoundaryCondition
 {

src/scheme/ExternalBoundaryConditionDescriptor.hpp

+#ifndef EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <mesh/IBoundaryDescriptor.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <utils/Socket.hpp>
+
+#include <memory>
+
+class ExternalBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "external(" << m_name << ',' << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  const std::string_view m_name;
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const std::shared_ptr<const Socket> m_socket;
+
+ public:
+  std::string_view
+  name() const
+  {
+    return m_name;
+  }
+
+  const std::shared_ptr<const Socket>&
+  socket() const
+  {
+    return m_socket;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const final
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::external;
+  }
+
+  ExternalBoundaryConditionDescriptor(const std::string_view name,
+                                      std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                      const std::shared_ptr<const Socket>& socket)
+    : m_name{name}, m_boundary_descriptor(boundary_descriptor), m_socket{socket}
+  {}
+
+  ExternalBoundaryConditionDescriptor(const ExternalBoundaryConditionDescriptor&) = delete;
+  ExternalBoundaryConditionDescriptor(ExternalBoundaryConditionDescriptor&&)      = delete;
+
+  ~ExternalBoundaryConditionDescriptor() = default;
+};
+
+#endif   // EXTERNAL_BOUNDARY_CONDITION_DESCRIPTOR_HPP

src/scheme/IBoundaryConditionDescriptor.hpp

@@ -12,6 +12,7 @@ class IBoundaryConditionDescriptor
   {
     axis,
     dirichlet,
+    external,
     fourier,
     fixed,
     free,