Add velocity boundary conditions for acoustic solver (Glace)

Still requires design improvements and to check data validity

src/language/modules/SchemeModule.cpp

@@ -10,6 +10,7 @@
 #include <scheme/NumberedBoundaryDescriptor.hpp>
 #include <scheme/PressureBoundaryConditionDescriptor.hpp>
 #include <scheme/SymmetryBoundaryConditionDescriptor.hpp>
+#include <scheme/VelocityBoundaryConditionDescriptor.hpp>
 
 #include <memory>
 
@@ -140,6 +141,50 @@ struct GlaceScheme
           }
           break;
         }
+        case IBoundaryConditionDescriptor::Type::velocity: {
+          const VelocityBoundaryConditionDescriptor& velocity_bc_descriptor =
+            dynamic_cast<const VelocityBoundaryConditionDescriptor&>(*bc_descriptor);
+          for (size_t i_ref_face_list = 0;
+               i_ref_face_list < m_mesh->connectivity().template numberOfRefItemList<FaceType>(); ++i_ref_face_list) {
+            const auto& ref_face_list = m_mesh->connectivity().template refItemList<FaceType>(i_ref_face_list);
+            const RefId& ref          = ref_face_list.refId();
+            if (ref == velocity_bc_descriptor.boundaryDescriptor()) {
+              const FunctionSymbolId velocity_id = velocity_bc_descriptor.functionSymbolId();
+
+              if constexpr (Dimension == 1) {
+                const auto& node_list = ref_face_list.list();
+
+                Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+                  TinyVector<Dimension>)>::template interpolate<FaceType>(velocity_id, m_mesh->xr(), node_list);
+
+                std::shared_ptr bc =
+                  std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
+                bc_list.push_back(BoundaryConditionHandler(bc));
+              } else {
+                const auto& face_list           = ref_face_list.list();
+                const auto& face_to_node_matrix = m_mesh->connectivity().faceToNodeMatrix();
+                std::set<NodeId> node_set;
+                for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+                  FaceId face_id         = face_list[i_face];
+                  const auto& face_nodes = face_to_node_matrix[face_id];
+                  for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+                    node_set.insert(face_nodes[i_node]);
+                  }
+                }
+
+                Array<NodeId> node_list = convert_to_array(node_set);
+
+                Array<const TinyVector<Dimension>> value_list = InterpolateItemValue<TinyVector<Dimension>(
+                  TinyVector<Dimension>)>::template interpolate<ItemType::node>(velocity_id, m_mesh->xr(), node_list);
+
+                std::shared_ptr bc =
+                  std::make_shared<VelocityBoundaryCondition<MeshType::Dimension>>(node_list, value_list);
+                bc_list.push_back(BoundaryConditionHandler(bc));
+              }
+            }
+          }
+          break;
+        }
         case IBoundaryConditionDescriptor::Type::pressure: {
           const PressureBoundaryConditionDescriptor& pressure_bc_descriptor =
             dynamic_cast<const PressureBoundaryConditionDescriptor&>(*bc_descriptor);
@@ -172,7 +217,9 @@ struct GlaceScheme
           break;
         }
         default: {
-          throw UnexpectedError("Unknown BCDescription\n");
+          std::ostringstream error_msg;
+          error_msg << *bc_descriptor << " is an invalid boundary condition for acoustic solver";
+          throw NormalError(error_msg.str());
         }
         }
       }
@@ -325,6 +372,19 @@ SchemeModule::SchemeModule()
 
                               ));
 
+  this->_addBuiltinFunction("velocity",
+                            std::make_shared<BuiltinFunctionEmbedder<std::shared_ptr<
+                              const IBoundaryConditionDescriptor>(std::shared_ptr<const IBoundaryDescriptor>,
+                                                                  const FunctionSymbolId&)>>(
+
+                              [](std::shared_ptr<const IBoundaryDescriptor> boundary,
+                                 const FunctionSymbolId& velocity_id)
+                                -> std::shared_ptr<const IBoundaryConditionDescriptor> {
+                                return std::make_shared<VelocityBoundaryConditionDescriptor>(boundary, velocity_id);
+                              }
+
+                              ));
+
   this->_addBuiltinFunction("glace",
                             std::make_shared<BuiltinFunctionEmbedder<
                               void(std::shared_ptr<const IMesh>,

src/scheme/AcousticSolver.hpp

@@ -123,11 +123,22 @@ class AcousticSolver
   {
     for (const auto& handler : m_boundary_condition_list) {
       switch (handler.boundaryCondition().type()) {
-      case BoundaryCondition::normal_velocity: {
-        throw NotImplementedError("normal_velocity BC");
-      }
       case BoundaryCondition::velocity: {
-        throw NotImplementedError("velocity BC");
+        const VelocityBoundaryCondition<Dimension>& velocity_bc =
+          dynamic_cast<const VelocityBoundaryCondition<Dimension>&>(handler.boundaryCondition());
+
+        const auto& node_list  = velocity_bc.nodeList();
+        const auto& value_list = velocity_bc.valueList();
+
+        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];
+
+            m_Ar[node_id] = identity;
+            m_br[node_id] = value;
+          });
+        break;
       }
       case BoundaryCondition::pressure: {
         const PressureBoundaryCondition<Dimension>& pressure_bc =
@@ -142,7 +153,8 @@ class AcousticSolver
 
           const auto& node_list  = pressure_bc.faceList();
           const auto& value_list = pressure_bc.valueList();
-          for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+          parallel_for(
+            node_list.size(), PUGS_LAMBDA(size_t i_node) {
               const NodeId node_id       = node_list[i_node];
               const auto& node_cell_list = node_to_cell_matrix[node_id];
               Assert(node_cell_list.size() == 1);
@@ -151,7 +163,7 @@ class AcousticSolver
               size_t node_local_number_in_cell = node_local_numbers_in_their_cells(node_id, 0);
 
               m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
-          }
+            });
         } else {
           const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
 

src/scheme/BoundaryCondition.hpp

@@ -15,7 +15,6 @@ class BoundaryCondition   // clazy:exclude=copyable-polymorphic
   enum Type
   {
     velocity,
-    normal_velocity,
     pressure,
     symmetry
   };
@@ -58,7 +57,7 @@ class PressureBoundaryCondition final : public BoundaryCondition   // clazy:excl
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const Array<const FaceId>& face_list, Array<const double> value_list)
+  PressureBoundaryCondition(const Array<const FaceId>& face_list, const Array<const double>& value_list)
     : BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, m_face_list{face_list}
   {}
 
@@ -85,7 +84,7 @@ class PressureBoundaryCondition<1> final : public BoundaryCondition   // clazy:e
     return m_value_list;
   }
 
-  PressureBoundaryCondition(const Array<const NodeId>& face_list, Array<const double> value_list)
+  PressureBoundaryCondition(const Array<const NodeId>& face_list, const Array<const double>& value_list)
     : BoundaryCondition{BoundaryCondition::pressure}, m_value_list{value_list}, m_face_list{face_list}
   {}
 
@@ -129,6 +128,33 @@ class SymmetryBoundaryCondition : public BoundaryCondition
   ~SymmetryBoundaryCondition() = default;
 };
 
+template <size_t Dimension>
+class VelocityBoundaryCondition final : public BoundaryCondition   // clazy:exclude=copyable-polymorphic
+{
+ private:
+  const Array<const TinyVector<Dimension>> m_value_list;
+  const Array<const NodeId> m_node_list;
+
+ public:
+  const Array<const NodeId>&
+  nodeList() const
+  {
+    return m_node_list;
+  }
+
+  const Array<const TinyVector<Dimension>>&
+  valueList() const
+  {
+    return m_value_list;
+  }
+
+  VelocityBoundaryCondition(const Array<const NodeId>& node_list, const Array<const TinyVector<Dimension>>& value_list)
+    : BoundaryCondition{BoundaryCondition::velocity}, m_value_list{value_list}, m_node_list{node_list}
+  {}
+
+  ~VelocityBoundaryCondition() = default;
+};
+
 class BoundaryConditionHandler
 {
  private:

src/scheme/IBoundaryConditionDescriptor.hpp

@@ -8,6 +8,7 @@ class IBoundaryConditionDescriptor
  public:
   enum class Type
   {
+    velocity,
     pressure,
     symmetry
   };

src/scheme/VelocityBoundaryConditionDescriptor.hpp

+#ifndef VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+#define VELOCITY_BOUNDARY_CONDITION_DESCRIPTOR_HPP
+
+#include <language/utils/FunctionSymbolId.hpp>
+#include <scheme/IBoundaryConditionDescriptor.hpp>
+#include <scheme/IBoundaryDescriptor.hpp>
+
+#include <memory>
+
+class VelocityBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
+{
+ private:
+  std::ostream&
+  _write(std::ostream& os) const final
+  {
+    os << "velocity(" << *m_boundary_descriptor << ")";
+    return os;
+  }
+
+  std::shared_ptr<const IBoundaryDescriptor> m_boundary_descriptor;
+  const FunctionSymbolId m_function_symbol_id;
+
+ public:
+  FunctionSymbolId
+  functionSymbolId() const
+  {
+    return m_function_symbol_id;
+  }
+
+  const IBoundaryDescriptor&
+  boundaryDescriptor() const
+  {
+    return *m_boundary_descriptor;
+  }
+
+  Type
+  type() const final
+  {
+    return Type::velocity;
+  }
+
+  VelocityBoundaryConditionDescriptor(std::shared_ptr<const IBoundaryDescriptor> boundary_descriptor,
+                                      const FunctionSymbolId& function_symbol_id)
+    : m_boundary_descriptor(boundary_descriptor), m_function_symbol_id{function_symbol_id}
+  {
+    ;
+  }
+
+  VelocityBoundaryConditionDescriptor(const VelocityBoundaryConditionDescriptor&) = delete;
+  VelocityBoundaryConditionDescriptor(VelocityBoundaryConditionDescriptor&&)      = delete;
+
+  ~VelocityBoundaryConditionDescriptor() = default;
+};
+
+#endif   // VELOCITY!_BOUNDARY_CONDITION_DESCRIPTOR_HPP