diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp index ed085fb72e3d2766a5b88105e5f32dd0bbba19b0..9f57841e16d414e2df6792705d521222250cbba0 100644 --- a/src/mesh/Connectivity.cpp +++ b/src/mesh/Connectivity.cpp @@ -150,6 +150,8 @@ void Connectivity<3>::_computeFaceCellConnectivities() m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_face_matrix, m_face_to_cell_matrix); + m_face_to_cell_local_face_matrix = CellValuePerFace<unsigned short>(*this); + m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_face_matrix, m_face_to_cell_matrix, m_face_to_cell_local_face_matrix); diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp index dd765a17314cc050d689d4c8bbb4bffba3c54c65..81dda065266c646bbd7430707f6425b1e3bb11ac 100644 --- a/src/mesh/Connectivity.hpp +++ b/src/mesh/Connectivity.hpp @@ -8,6 +8,7 @@ #include <ConnectivityMatrix.hpp> #include <ConnectivityComputer.hpp> +#include <SubItemValuePerItem.hpp> #include <vector> #include <unordered_map> @@ -21,6 +22,8 @@ #include <tuple> #include <algorithm> +#include <IConnectivity.hpp> + template <size_t Dimension> class Connectivity; @@ -218,7 +221,8 @@ class ConnectivityFace<3> }; template <size_t Dimension> -class Connectivity +class Connectivity final + : public IConnectivity { public: static constexpr size_t dimension = Dimension; @@ -231,16 +235,19 @@ public: ConnectivityMatrixShort m_cell_to_face_is_reversed_matrix; ConnectivityMatrix m_face_to_cell_matrix; - ConnectivityMatrixShort m_face_to_cell_local_face_matrix; + CellValuePerFace<unsigned short> m_face_to_cell_local_face_matrix; ConnectivityMatrix m_face_to_node_matrix; ConnectivityMatrix m_node_to_cell_matrix; - ConnectivityMatrixShort m_node_to_cell_local_node_matrix; + CellValuePerNode<unsigned short> m_node_to_cell_local_node_matrix; template <TypeOfItem SubItemType, TypeOfItem ItemType> const ConnectivityMatrix& itemToItemMatrix() const = delete; + const ConnectivityMatrix& itemToItemMatrix(const TypeOfItem& item_type_0, + const TypeOfItem& item_type_1) const final; + private: ConnectivityComputer m_connectivity_computer; @@ -344,6 +351,8 @@ private: m_connectivity_computer.computeInverseConnectivityMatrix(m_cell_to_node_matrix, m_node_to_cell_matrix); + m_node_to_cell_local_node_matrix = CellValuePerNode<unsigned short>(*this); + m_connectivity_computer.computeLocalChildItemNumberInItem(m_cell_to_node_matrix, m_node_to_cell_matrix, m_node_to_cell_local_node_matrix); @@ -361,6 +370,15 @@ private: using Connectivity3D = Connectivity<3>; +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<3>::itemToItemMatrix<TypeOfItem::cell, + TypeOfItem::face>() const +{ + return m_cell_to_face_matrix; +} + template <> template <> inline const ConnectivityMatrix& @@ -370,8 +388,36 @@ Connectivity<3>::itemToItemMatrix<TypeOfItem::cell, return m_cell_to_node_matrix; } +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<3>::itemToItemMatrix<TypeOfItem::face, + TypeOfItem::cell>() const +{ + return m_face_to_cell_matrix; +} + +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<3>::itemToItemMatrix<TypeOfItem::node, + TypeOfItem::cell>() const +{ + return m_node_to_cell_matrix; +} + + using Connectivity2D = Connectivity<2>; +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<2>::itemToItemMatrix<TypeOfItem::cell, + TypeOfItem::face>() const +{ + return m_cell_to_face_matrix; +} + template <> template <> inline const ConnectivityMatrix& @@ -381,6 +427,24 @@ Connectivity<2>::itemToItemMatrix<TypeOfItem::cell, return m_cell_to_node_matrix; } +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<2>::itemToItemMatrix<TypeOfItem::face, + TypeOfItem::cell>() const +{ + return m_face_to_cell_matrix; +} + +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<2>::itemToItemMatrix<TypeOfItem::node, + TypeOfItem::cell>() const +{ + return m_node_to_cell_matrix; +} + using Connectivity1D = Connectivity<1>; template <> @@ -392,4 +456,82 @@ Connectivity<1>::itemToItemMatrix<TypeOfItem::cell, return m_cell_to_node_matrix; } +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<1>::itemToItemMatrix<TypeOfItem::cell, + TypeOfItem::face>() const +{ +#warning in 1d, faces and node are the same + return m_cell_to_face_matrix; +} + +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<1>::itemToItemMatrix<TypeOfItem::face, + TypeOfItem::cell>() const +{ +#warning in 1d, faces and node are the same + return m_face_to_cell_matrix; +} + +template <> +template <> +inline const ConnectivityMatrix& +Connectivity<1>::itemToItemMatrix<TypeOfItem::node, + TypeOfItem::cell>() const +{ +#warning in 1d, faces and node are the same + return m_node_to_cell_matrix; +} + +template <size_t Dimension> +const ConnectivityMatrix& +Connectivity<Dimension>:: +itemToItemMatrix(const TypeOfItem& item_type_0, + const TypeOfItem& item_type_1) const +{ + switch (item_type_0) { + case TypeOfItem::cell: { + switch (item_type_1) { + case TypeOfItem::node: { + return itemToItemMatrix<TypeOfItem::cell, TypeOfItem::node>(); + } + default: { + std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n"; + std::exit(1); + } + } + } + case TypeOfItem::face: { + switch (item_type_1) { + case TypeOfItem::cell: { + return itemToItemMatrix<TypeOfItem::face, TypeOfItem::cell>(); + } + default: { + std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n"; + std::exit(1); + } + } + } + case TypeOfItem::node: { + switch (item_type_1) { + case TypeOfItem::cell: { + return itemToItemMatrix<TypeOfItem::node, TypeOfItem::cell>(); + } + default: { + std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_1) << "\n"; + std::exit(1); + } + } + } + default: { + std::cerr << __FILE__ << ":" << __LINE__ << ": NIY " << int(item_type_0) << "\n"; + std::exit(1); + } + } +} + + #endif // CONNECTIVITY_HPP diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp index b950261546d6f24c4a8612de54f8950ff9417a19..ea6900fa1546591b77e2f39c1211aba9f467f344 100644 --- a/src/mesh/ConnectivityComputer.cpp +++ b/src/mesh/ConnectivityComputer.cpp @@ -36,23 +36,43 @@ computeInverseConnectivityMatrix(const ConnectivityMatrix& item_to_child_item_ma void ConnectivityComputer:: computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix, const ConnectivityMatrix& child_item_to_items_matrix, - ConnectivityMatrixShort& child_item_number_in_item_matrix) const + CellValuePerNode<unsigned short>& child_item_number_in_item_matrix) const { - std::vector<std::vector<unsigned int>> child_item_number_in_item_vector(child_item_to_items_matrix.numRows()); for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) { const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r); - child_item_number_in_item_vector[r].resize(child_item_to_items.length); for (unsigned short J=0; J<child_item_to_items.length; ++J) { const unsigned int j = child_item_to_items(J); const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j); for (unsigned int R=0; R<item_to_child_items.length; ++R) { if (item_to_child_items(R) == r) { - child_item_number_in_item_vector[r][J] = R; + child_item_number_in_item_matrix(r,J) = R; + break; + } + } + } + } +} + + + +void ConnectivityComputer:: +computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix, + const ConnectivityMatrix& child_item_to_items_matrix, + CellValuePerFace<unsigned short>& child_item_number_in_item_matrix) const +{ + for (unsigned int r=0; r<child_item_to_items_matrix.numRows(); ++r) { + const auto& child_item_to_items = child_item_to_items_matrix.rowConst(r); + for (unsigned short J=0; J<child_item_to_items.length; ++J) { + const unsigned int j = child_item_to_items(J); + const auto& item_to_child_items = item_to_child_items_matrix.rowConst(j); + + for (unsigned int R=0; R<item_to_child_items.length; ++R) { + if (item_to_child_items(R) == r) { + child_item_number_in_item_matrix(r,J) = R; break; } } } } - child_item_number_in_item_matrix = child_item_number_in_item_vector; } diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp index 1b8ee0180affa5007feb7d7482304e1a54ae3a93..4fb47cc206edcc34e5186e9e37dac03b15741cd7 100644 --- a/src/mesh/ConnectivityComputer.hpp +++ b/src/mesh/ConnectivityComputer.hpp @@ -2,6 +2,7 @@ #define CONNECTIVITY_COMPUTER_HPP #include <ConnectivityMatrix.hpp> +#include <SubItemValuePerItem.hpp> struct ConnectivityComputer { @@ -10,7 +11,11 @@ struct ConnectivityComputer void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix, const ConnectivityMatrix& child_item_to_items_matrix, - ConnectivityMatrixShort& child_item_number_in_item_matrix) const; + CellValuePerNode<unsigned short>& child_item_number_in_item_matrix) const; + + void computeLocalChildItemNumberInItem(const ConnectivityMatrix& item_to_child_items_matrix, + const ConnectivityMatrix& child_item_to_items_matrix, + CellValuePerFace<unsigned short>& child_item_number_in_item_matrix) const; }; #endif // CONNECTIVITY_COMPUTER_HPP diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp new file mode 100644 index 0000000000000000000000000000000000000000..6dc6ca4a2586f4fce344614f39c5bc82b69dd64c --- /dev/null +++ b/src/mesh/IConnectivity.hpp @@ -0,0 +1,18 @@ +#ifndef ICONNECTIVITY_HPP +#define ICONNECTIVITY_HPP + +#include <TypeOfItem.hpp> +#include <ConnectivityMatrix.hpp> + +struct IConnectivity +{ + virtual const ConnectivityMatrix& + itemToItemMatrix(const TypeOfItem& item_type_0, + const TypeOfItem& item_type_1) const = 0; + + IConnectivity() = default; + IConnectivity(const IConnectivity&) = delete; + virtual ~IConnectivity() = default; +}; + +#endif // ICONNECTIVITY_HPP diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp index a61c863cfd5687db91aefb8f9f3f012c75f4609f..384e91c7fc92e1f79c9ad51f774332cdabc17768 100644 --- a/src/mesh/SubItemValuePerItem.hpp +++ b/src/mesh/SubItemValuePerItem.hpp @@ -3,8 +3,10 @@ #include <Kokkos_StaticCrsGraph.hpp> #include <TypeOfItem.hpp> +#include <PastisAssert.hpp> #include <ConnectivityMatrix.hpp> +#include <IConnectivity.hpp> template <typename DataType, TypeOfItem SubItemType, @@ -158,11 +160,10 @@ class SubItemValuePerItem SubItemValuePerItem() = default; - template <typename ConnectivityType> - SubItemValuePerItem(const ConnectivityType& connectivity) + SubItemValuePerItem(const IConnectivity& connectivity) { ConnectivityMatrix connectivity_matrix - = connectivity.template itemToItemMatrix<ItemType,SubItemType>(); + = connectivity.itemToItemMatrix(ItemType,SubItemType); m_host_row_map = connectivity_matrix.rowsMap(); m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries()); } @@ -173,4 +174,10 @@ class SubItemValuePerItem template <typename DataType> using NodeValuePerCell = SubItemValuePerItem<DataType, TypeOfItem::node, TypeOfItem::cell>; +template <typename DataType> +using CellValuePerNode = SubItemValuePerItem<DataType, TypeOfItem::cell, TypeOfItem::node>; + +template <typename DataType> +using CellValuePerFace = SubItemValuePerItem<DataType, TypeOfItem::cell, TypeOfItem::face>; + #endif // SUBITEM_VALUE_PER_ITEM_HPP diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp index b00571b597ffb9401bcd4098f846492c2c65571c..99e1743140e5bd11fe334f88e14c99b7bebe58d4 100644 --- a/src/scheme/AcousticSolver.hpp +++ b/src/scheme/AcousticSolver.hpp @@ -101,10 +101,10 @@ class AcousticSolver Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) { Rdd sum = zero; const auto& node_to_cell = m_connectivity.m_node_to_cell_matrix.rowConst(r); - const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.rowConst(r); + const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.itemValues(r); for (size_t j=0; j<node_to_cell.length; ++j) { const unsigned int J = node_to_cell(j); - const unsigned int R = node_to_cell_local_node(j); + const unsigned int R = node_to_cell_local_node[j]; sum += Ajr(J,R); } m_Ar(r) = sum; @@ -123,10 +123,10 @@ class AcousticSolver Rd& br = m_br(r); br = zero; const auto& node_to_cell = m_connectivity.m_node_to_cell_matrix.rowConst(r); - const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.rowConst(r); + const auto& node_to_cell_local_node = m_connectivity.m_node_to_cell_local_node_matrix.itemValues(r); for (size_t j=0; j<node_to_cell.length; ++j) { const unsigned int J = node_to_cell(j); - const unsigned int R = node_to_cell_local_node(j); + const unsigned int R = node_to_cell_local_node[j]; br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R); } });