Add pressure boundary conditions for acoustic solver

src/language/modules/SchemeModule.cpp

@@ -57,6 +57,39 @@ class InterpolateItemValue<OutputType(InputType)> : public PugsFunctionAdapter<O
 
     return value;
   }
+
+  template <ItemType item_type>
+  static inline Array<OutputType>
+  interpolate(const FunctionSymbolId& function_symbol_id,
+              const ItemValue<const InputType, item_type>& position,
+              const Array<const ItemIdT<item_type>>& list_of_items)
+  {
+    auto& expression    = Adapter::getFunctionExpression(function_symbol_id);
+    auto convert_result = Adapter::getResultConverter(expression.m_data_type);
+
+    Array<ExecutionPolicy> context_list = Adapter::getContextList(expression);
+
+    using execution_space = typename Kokkos::DefaultExecutionSpace::execution_space;
+    Kokkos::Experimental::UniqueToken<execution_space, Kokkos::Experimental::UniqueTokenScope::Global> tokens;
+
+    Array<OutputType> value{list_of_items.size()};
+    using ItemId = ItemIdT<item_type>;
+
+    parallel_for(list_of_items.size(), [=, &expression, &tokens](size_t i_item) {
+      ItemId item_id  = list_of_items[i_item];
+      const int32_t t = tokens.acquire();
+
+      auto& execution_policy = context_list[t];
+
+      Adapter::convertArgs(execution_policy.currentContext(), position[item_id]);
+      auto result   = expression.execute(execution_policy);
+      value[i_item] = convert_result(std::move(result));
+
+      tokens.release(t);
+    });
+
+    return value;
+  }
 };
 
 template <size_t Dimension>
@@ -116,8 +149,20 @@ struct GlaceScheme
             const RefId& ref          = ref_face_list.refId();
             if (ref == pressure_bc_descriptor.boundaryDescriptor()) {
               const auto& face_list = ref_face_list.list();
-              Array<double> face_values{face_list.size()};
-              face_values.fill(1);
+
+              const FunctionSymbolId pressure_id = pressure_bc_descriptor.functionSymbolId();
+
+              Array<const double> face_values = [&] {
+                if constexpr (Dimension == 1) {
+                  return InterpolateItemValue<double(
+                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, m_mesh->xr(), face_list);
+                } else {
+                  MeshDataType& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+
+                  return InterpolateItemValue<double(
+                    TinyVector<Dimension>)>::template interpolate<FaceType>(pressure_id, mesh_data.xl(), face_list);
+                }
+              }();
 
               std::shared_ptr bc =
                 std::make_shared<PressureBoundaryCondition<MeshType::Dimension>>(face_list, face_values);

src/mesh/MeshData.hpp

@@ -35,13 +35,37 @@ class MeshData : public IMeshData
   std::shared_ptr<NodeValuePerCell<const Rd>> m_Cjr;
   std::shared_ptr<NodeValuePerCell<const double>> m_ljr;
   std::shared_ptr<NodeValuePerCell<const Rd>> m_njr;
+  FaceValue<const Rd> m_xl;
   std::shared_ptr<CellValue<const Rd>> m_xj;
   std::shared_ptr<CellValue<const double>> m_Vj;
   std::shared_ptr<NodeValuePerFace<const Rd>> m_Nlr;
 
   PUGS_INLINE
   void
-  _computeIsobarycenter()
+  _computeFaceIsobarycenter()
+  {   // Computes vertices isobarycenter
+    if constexpr (Dimension == 1) {
+      static_assert(Dimension != 1, "xl does not make sense in 1d");
+    } else {
+      const NodeValue<const Rd>& xr = m_mesh.xr();
+
+      const auto& face_to_node_matrix = m_mesh.connectivity().faceToNodeMatrix();
+      FaceValue<Rd> xl(m_mesh.connectivity());
+      parallel_for(
+        m_mesh.numberOfFaces(), PUGS_LAMBDA(FaceId face_id) {
+          Rd X                   = zero;
+          const auto& face_nodes = face_to_node_matrix[face_id];
+          for (size_t R = 0; R < face_nodes.size(); ++R) {
+            X += xr[face_nodes[R]];
+          }
+          xl[face_id] = 1. / face_nodes.size() * X;
+        });
+      m_xl = xl;
+    }
+  }
+
+  void
+  _computeCellIsobarycenter()
   {   // Computes vertices isobarycenter
     if constexpr (Dimension == 1) {
       const NodeValue<const Rd>& xr = m_mesh.xr();
@@ -337,12 +361,22 @@ class MeshData : public IMeshData
     return *m_njr;
   }
 
+  PUGS_INLINE
+  FaceValue<const Rd>
+  xl()
+  {
+    if (not m_xl.isBuilt()) {
+      this->_computeFaceIsobarycenter();
+    }
+    return m_xl;
+  }
+
   PUGS_INLINE
   CellValue<const Rd>
   xj()
   {
     if (not m_xj) {
-      this->_computeIsobarycenter();
+      this->_computeCellIsobarycenter();
     }
     return *m_xj;
   }

src/scheme/AcousticSolver.hpp

@@ -132,8 +132,9 @@ class AcousticSolver
       case BoundaryCondition::pressure: {
         const PressureBoundaryCondition<Dimension>& pressure_bc =
           dynamic_cast<const PressureBoundaryCondition<Dimension>&>(handler.boundaryCondition());
-        if constexpr (Dimension == 1) {
+
         MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
+        if constexpr (Dimension == 1) {
           const NodeValuePerCell<const Rd> Cjr = mesh_data.Cjr();
 
           const auto& node_to_cell_matrix               = m_connectivity.nodeToCellMatrix();
@@ -152,7 +153,32 @@ class AcousticSolver
             m_br[node_id] -= value_list[i_node] * Cjr(node_cell_id, node_local_number_in_cell);
           }
         } else {
-          throw NotImplementedError("pressure bc in dimension>1");
+          const NodeValuePerFace<const Rd> Nlr = mesh_data.Nlr();
+
+          const auto& face_to_cell_matrix               = m_connectivity.faceToCellMatrix();
+          const auto& face_to_node_matrix               = m_connectivity.faceToNodeMatrix();
+          const auto& face_local_numbers_in_their_cells = m_connectivity.faceLocalNumbersInTheirCells();
+          const auto& face_cell_is_reversed             = m_connectivity.cellFaceIsReversed();
+
+          const auto& face_list  = pressure_bc.faceList();
+          const auto& value_list = pressure_bc.valueList();
+          for (size_t i_face = 0; i_face < face_list.size(); ++i_face) {
+            const FaceId face_id       = face_list[i_face];
+            const auto& face_cell_list = face_to_cell_matrix[face_id];
+            Assert(face_cell_list.size() == 1);
+
+            CellId face_cell_id              = face_cell_list[0];
+            size_t face_local_number_in_cell = face_local_numbers_in_their_cells(face_id, 0);
+
+            const double sign = face_cell_is_reversed(face_cell_id, face_local_number_in_cell) ? -1 : 1;
+
+            const auto& face_nodes = face_to_node_matrix[face_id];
+
+            for (size_t i_node = 0; i_node < face_nodes.size(); ++i_node) {
+              NodeId node_id = face_nodes[i_node];
+              m_br[node_id] -= sign * value_list[i_face] * Nlr(face_id, i_node);
+            }
+          }
         }
         break;
       }

src/scheme/PressureBoundaryConditionDescriptor.hpp

@@ -21,6 +21,12 @@ class PressureBoundaryConditionDescriptor : public IBoundaryConditionDescriptor
   const FunctionSymbolId m_function_symbol_id;
 
  public:
+  FunctionSymbolId
+  functionSymbolId() const
+  {
+    return m_function_symbol_id;
+  }
+
   const IBoundaryDescriptor&
   boundaryDescriptor() const
   {