diff --git a/src/language/modules/SchemeModule.cpp b/src/language/modules/SchemeModule.cpp
index 8ff18d22f45c273e38ed907d5ff23fa9809ae7e8..7903e6ce4b18fd987de1820fa9b164ddc99d352e 100644
--- a/src/language/modules/SchemeModule.cpp
+++ b/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);
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 86e78f5f7fde25656882535258f6fb3779855d01..d3229925150a709056547abcd6550670c2e2e6d0 100644
--- a/src/mesh/MeshData.hpp
+++ b/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;
}
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 622b998931ff08f9db141469da71b6f1aaf0105b..75b8eec54f8da6b9f97d30789b98044436354287 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/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());
+
+ MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
if constexpr (Dimension == 1) {
- MeshData<Dimension>& mesh_data = MeshDataManager::instance().getMeshData(*m_mesh);
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;
}
diff --git a/src/scheme/PressureBoundaryConditionDescriptor.hpp b/src/scheme/PressureBoundaryConditionDescriptor.hpp
index 28cdb204e1db9a87c52b84aedc7c1ee9b0b503e6..8a3d39f5138e05fa41e856ebd676a0cc0c729767 100644
--- a/src/scheme/PressureBoundaryConditionDescriptor.hpp
+++ b/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
{