From 1a6edd28a9f95556b376bfbd0b336660d84847fb Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 11 Sep 2018 18:44:18 +0200
Subject: [PATCH] Improve access to item->item connectivity
Returned type is now ItemId
---
src/mesh/Connectivity.cpp | 4 +-
src/mesh/Connectivity.hpp | 115 ++++++++++++++++++++++++++----
src/mesh/ConnectivityComputer.cpp | 6 +-
src/mesh/IConnectivity.hpp | 15 +++-
src/mesh/ItemToItemMatrix.hpp | 77 ++++++++++++++++++++
src/mesh/MeshData.hpp | 84 ++++++++++------------
src/mesh/MeshNodeBoundary.hpp | 57 +++++++--------
src/mesh/SubItemValuePerItem.hpp | 2 +-
src/output/VTKWriter.hpp | 23 +++---
src/scheme/AcousticSolver.hpp | 43 +++++------
10 files changed, 290 insertions(+), 136 deletions(-)
create mode 100644 src/mesh/ItemToItemMatrix.hpp
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index c0bfbc541..f84f3835a 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -7,7 +7,7 @@ void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell, ItemType::node);
+ = this->_getMatrix(ItemType::cell, ItemType::node);
CellValue<unsigned short> cell_nb_faces(*this);
std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
@@ -157,7 +157,7 @@ template<>
void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
{
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell, ItemType::node);
+ = this->_getMatrix(ItemType::cell, ItemType::node);
// In 2D faces are simply define
using CellFaceId = std::pair<CellId, unsigned short>;
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index c2784d337..cfcb7f679 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -10,6 +10,8 @@
#include <IConnectivity.hpp>
#include <ConnectivityMatrix.hpp>
+#include <ItemToItemMatrix.hpp>
+
#include <SubItemValuePerItem.hpp>
#include <ConnectivityComputer.hpp>
@@ -223,6 +225,7 @@ class ConnectivityFace<3>
~ConnectivityFace() = default;
};
+
template <size_t Dimension>
class Connectivity final
: public IConnectivity
@@ -282,6 +285,23 @@ class Connectivity final
return sub_item_value_per_item;
}
+ friend class ConnectivityComputer;
+
+ KOKKOS_INLINE_FUNCTION
+ const ConnectivityMatrix& _getMatrix(const ItemType& item_type_0,
+ const ItemType& item_type_1) const
+ {
+ const ConnectivityMatrix& connectivity_matrix
+ = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
+ if (not connectivity_matrix.isBuilt()) {
+ const_cast<ConnectivityMatrix&>(connectivity_matrix)
+ = m_connectivity_computer
+ . computeConnectivityMatrix(*this, item_type_0, item_type_1);
+
+ }
+ return connectivity_matrix;
+ }
+
public:
KOKKOS_INLINE_FUNCTION
@@ -292,22 +312,89 @@ class Connectivity final
= m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
return connectivity_matrix.isBuilt();
}
+ template <ItemType source_item_type,
+ ItemType target_item_type>
+ KOKKOS_INLINE_FUNCTION
+ const auto getItemToItemMatrix() const
+ {
+ return ItemToItemMatrix<source_item_type,
+ target_item_type>(_getMatrix(source_item_type,
+ target_item_type));
+ }
KOKKOS_INLINE_FUNCTION
- const ConnectivityMatrix& getMatrix(const ItemType& item_type_0,
- const ItemType& item_type_1) const
+ const auto cellToFaceMatrix() const
{
- const ConnectivityMatrix& connectivity_matrix
- = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
- if (not connectivity_matrix.isBuilt()) {
- const_cast<ConnectivityMatrix&>(connectivity_matrix)
- = m_connectivity_computer
- . computeConnectivityMatrix(*this, item_type_0, item_type_1);
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::face>();
+ }
- }
- return connectivity_matrix;
+ KOKKOS_INLINE_FUNCTION
+ const auto cellToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::edge>();
}
+ KOKKOS_INLINE_FUNCTION
+ const auto cellToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::cell, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::edge>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto faceToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::face, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToFaceMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::face>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto edgeToNodeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::edge, ItemType::node>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToCellMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::cell>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToFaceMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::face>();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto nodeToEdgeMatrix() const
+ {
+ return this->template getItemToItemMatrix<ItemType::node, ItemType::edge>();
+ }
+
+
KOKKOS_INLINE_FUNCTION
const auto& cellFaceIsReversed() const
{
@@ -443,7 +530,7 @@ class Connectivity final
size_t numberOfNodes() const final
{
const auto& node_to_cell_matrix
- = this->getMatrix(ItemType::node,ItemType::cell);
+ = this->_getMatrix(ItemType::node,ItemType::cell);
return node_to_cell_matrix.numRows();
}
@@ -451,7 +538,7 @@ class Connectivity final
size_t numberOfEdges() const final
{
const auto& edge_to_node_matrix
- = this->getMatrix(ItemType::edge,ItemType::node);
+ = this->_getMatrix(ItemType::edge,ItemType::node);
return edge_to_node_matrix.numRows();
}
@@ -459,7 +546,7 @@ class Connectivity final
size_t numberOfFaces() const final
{
const auto& face_to_node_matrix
- = this->getMatrix(ItemType::face,ItemType::cell);
+ = this->_getMatrix(ItemType::face,ItemType::cell);
return face_to_node_matrix.numRows();
}
@@ -467,7 +554,7 @@ class Connectivity final
size_t numberOfCells() const final
{
const auto& cell_to_node_matrix
- = this->getMatrix(ItemType::cell,ItemType::node);
+ = this->_getMatrix(ItemType::cell,ItemType::node);
return cell_to_node_matrix.numRows();
}
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index 672502070..aa6f5c56a 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -14,7 +14,7 @@ computeConnectivityMatrix(const ConnectivityType& connectivity,
ConnectivityMatrix item_to_child_item_matrix;
if (connectivity.isConnectivityMatrixBuilt(child_item_type, item_type)) {
const ConnectivityMatrix& child_to_item_matrix
- = connectivity.getMatrix(child_item_type, item_type);
+ = connectivity._getMatrix(child_item_type, item_type);
std::cout << "computing connectivity "
<< itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
@@ -83,9 +83,9 @@ SubItemValuePerItem<const unsigned short, child_item_type, item_type>
ConnectivityComputer::computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const
{
const ConnectivityMatrix& child_item_to_items_matrix
- = connectivity.getMatrix(child_item_type, item_type);
+ = connectivity._getMatrix(child_item_type, item_type);
const ConnectivityMatrix& item_to_child_items_matrix
- = connectivity.getMatrix(item_type, child_item_type);
+ = connectivity._getMatrix(item_type, child_item_type);
SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
for (ItemIdT<item_type> r=0; r<item_to_child_items_matrix.numRows(); ++r) {
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
index 27ae9fa81..e27c2ef08 100644
--- a/src/mesh/IConnectivity.hpp
+++ b/src/mesh/IConnectivity.hpp
@@ -4,12 +4,21 @@
#include <ItemType.hpp>
#include <ConnectivityMatrix.hpp>
-struct IConnectivity
+class IConnectivity
{
+ protected:
+ template <typename DataType,
+ ItemType sub_item_type,
+ ItemType item_type,
+ typename Allowed>
+ friend
+ class SubItemValuePerItem;
+
virtual const ConnectivityMatrix&
- getMatrix(const ItemType& item_type_0,
- const ItemType& item_type_1) const = 0;
+ _getMatrix(const ItemType& item_type_0,
+ const ItemType& item_type_1) const = 0;
+ public:
virtual size_t numberOfNodes() const = 0;
virtual size_t numberOfEdges() const = 0;
virtual size_t numberOfFaces() const = 0;
diff --git a/src/mesh/ItemToItemMatrix.hpp b/src/mesh/ItemToItemMatrix.hpp
new file mode 100644
index 000000000..d03ee9e0c
--- /dev/null
+++ b/src/mesh/ItemToItemMatrix.hpp
@@ -0,0 +1,77 @@
+#ifndef ITEM_TO_ITEM_MATRIX_HPP
+#define ITEM_TO_ITEM_MATRIX_HPP
+
+#include <ItemType.hpp>
+#include <ItemId.hpp>
+
+#include <Kokkos_Core.hpp>
+#include <ConnectivityMatrix.hpp>
+
+template <ItemType SourceItemType,
+ ItemType TargetItemType>
+class ItemToItemMatrix
+{
+ public:
+ using SourceItemId = ItemIdT<SourceItemType>;
+ using TargetItemId = ItemIdT<TargetItemType>;
+
+ template <typename RowType>
+ class SubItemList
+ {
+ private:
+ const RowType m_row;
+
+ public:
+
+ KOKKOS_INLINE_FUNCTION
+ size_t size() const
+ {
+ return m_row.length;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ TargetItemId operator[](const size_t& j) const
+ {
+ return m_row(j);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ SubItemList(const RowType& row)
+ : m_row{row}
+ {
+ ;
+ }
+
+ SubItemList& operator=(const SubItemList&) = default;
+ SubItemList& operator=(SubItemList&&) = default;
+
+ SubItemList(const SubItemList&) = default;
+ SubItemList(SubItemList&&) = default;
+ ~SubItemList() = default;
+ };
+
+ private:
+ const ConnectivityMatrix& m_connectivity_matrix;
+
+ public:
+ const auto operator[](const SourceItemId& source_id) const
+ {
+ using RowType = decltype(m_connectivity_matrix.rowConst(source_id));
+ return SubItemList<RowType>(m_connectivity_matrix.rowConst(source_id));
+ }
+
+ ItemToItemMatrix(const ConnectivityMatrix& connectivity_matrix)
+ : m_connectivity_matrix{connectivity_matrix}
+ {
+ ;
+ }
+
+ ItemToItemMatrix& operator=(const ItemToItemMatrix&) = default;
+ ItemToItemMatrix& operator=(ItemToItemMatrix&&) = default;
+
+ ItemToItemMatrix(ItemToItemMatrix&&) = default;
+ ItemToItemMatrix(const ItemToItemMatrix&) = default;
+ ~ItemToItemMatrix() = default;
+};
+
+#endif // ITEM_TO_ITEM_MATRIX_HPP
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 2c8689926..1a07f1bcc 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -37,13 +37,12 @@ class MeshData
const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
CellValue<Rd> xj(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- xj[j] = 0.5*(xr[NodeId{cell_nodes(0)}]+xr[NodeId{cell_nodes(1)}]);
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ xj[j] = 0.5*(xr[cell_nodes[0]]+xr[cell_nodes[1]]);
});
m_xj = xj;
} else {
@@ -53,14 +52,13 @@ class MeshData
= m_mesh.connectivity().invCellNbNodes();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
CellValue<Rd> xj(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
Rd X = zero;
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (size_t R=0; R<cell_nodes.length; ++R) {
- X += xr[NodeId{cell_nodes(R)}];
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ X += xr[cell_nodes[R]];
}
xj[j] = inv_cell_nb_nodes[j]*X;
});
@@ -73,15 +71,15 @@ class MeshData
{
const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
+
CellValue<double> Vj(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
double sum_cjr_xr = 0;
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t R=0; R<cell_nodes.length; ++R) {
- sum_cjr_xr += (xr[NodeId{cell_nodes(R)}], m_Cjr(j,R));
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ sum_cjr_xr += (xr[cell_nodes[R]], m_Cjr(j,R));
}
Vj[j] = inv_dimension * sum_cjr_xr;
});
@@ -96,17 +94,16 @@ class MeshData
else if constexpr (dimension == 2) {
const NodeValue<const Rd>& xr = m_mesh.xr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
{
NodeValuePerCell<Rd> Cjr(m_mesh.connectivity());
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (size_t R=0; R<cell_nodes.length; ++R) {
- int Rp1 = (R+1)%cell_nodes.length;
- int Rm1 = (R+cell_nodes.length-1)%cell_nodes.length;
- Rd half_xrp_xrm = 0.5*(xr[NodeId{cell_nodes(Rp1)}]-xr[NodeId{cell_nodes(Rm1)}]);
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ int Rp1 = (R+1)%cell_nodes.size();
+ int Rm1 = (R+cell_nodes.size()-1)%cell_nodes.size();
+ Rd half_xrp_xrm = 0.5*(xr[cell_nodes[Rp1]]-xr[cell_nodes[Rm1]]);
Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
}
});
@@ -133,38 +130,35 @@ class MeshData
NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
const auto& face_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::face,
- ItemType::node);
+ = m_mesh.connectivity().faceToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const FaceId& l) {
- const auto& face_nodes = face_to_node_matrix.rowConst(l);
- const size_t nb_nodes = face_nodes.length;
+ const auto& face_nodes = face_to_node_matrix[l];
+ const size_t nb_nodes = face_nodes.size();
std::vector<Rd> dxr(nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
dxr[r]
- = xr[NodeId{face_nodes((r+1)%nb_nodes)}]
- - xr[NodeId{face_nodes((r+nb_nodes-1)%nb_nodes)}];
+ = xr[face_nodes[(r+1)%nb_nodes]]
+ - xr[face_nodes[(r+nb_nodes-1)%nb_nodes]];
}
const double inv_12_nb_nodes = 1./(12.*nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
Rd Nr = zero;
const Rd two_dxr = 2*dxr[r];
for (size_t s=0; s<nb_nodes; ++s) {
- Nr += crossProduct((two_dxr - dxr[s]), xr[NodeId{face_nodes(s)}]);
+ Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes[s]]);
}
Nr *= inv_12_nb_nodes;
- Nr -= (1./6.)*crossProduct(dxr[r], xr[NodeId{face_nodes(r)}]);
+ Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes[r]]);
Nlr(l,r) = Nr;
}
});
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
const auto& cell_to_face_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::face);
+ = m_mesh.connectivity().cellToFaceMatrix();
const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
@@ -175,20 +169,20 @@ class MeshData
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
- const auto& cell_faces = cell_to_face_matrix.rowConst(j);
+ const auto& cell_faces = cell_to_face_matrix[j];
const auto& face_is_reversed = cell_face_is_reversed.itemValues(j);
- for (size_t L=0; L<cell_faces.length; ++L) {
- const size_t l = cell_faces(L);
- const auto& face_nodes = face_to_node_matrix.rowConst(l);
+ for (size_t L=0; L<cell_faces.size(); ++L) {
+ const FaceId& l = cell_faces[L];
+ const auto& face_nodes = face_to_node_matrix[l];
#warning should this lambda be replaced by a precomputed correspondance?
std::function local_node_number_in_cell
- = [&](const size_t& node_number) {
- for (size_t i_node=0; i_node<cell_nodes.length; ++i_node) {
- if (node_number == cell_nodes(i_node)) {
+ = [&](const NodeId& node_number) {
+ for (size_t i_node=0; i_node<cell_nodes.size(); ++i_node) {
+ if (node_number == cell_nodes[i_node]) {
return i_node;
break;
}
@@ -197,13 +191,13 @@ class MeshData
};
if (face_is_reversed[L]) {
- for (size_t rl = 0; rl<face_nodes.length; ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes(rl));
+ for (size_t rl = 0; rl<face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
Cjr(j, R) -= Nlr(l,rl);
}
} else {
- for (size_t rl = 0; rl<face_nodes.length; ++rl) {
- const size_t R = local_node_number_in_cell(face_nodes(rl));
+ for (size_t rl = 0; rl<face_nodes.size(); ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes[rl]);
Cjr(j, R) += Nlr(l,rl);
}
}
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 8500e2e9f..fe84d7860 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -34,12 +34,12 @@ class MeshNodeBoundary
{
static_assert(dimension == MeshType::dimension);
const auto& face_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::face, ItemType::cell);
+ = mesh.connectivity().faceToCellMatrix();
const Array<const FaceId>& face_list = ref_face_list.faceList();
Kokkos::parallel_for(face_list.size(), KOKKOS_LAMBDA(const int& l){
- const auto& face_cells = face_to_cell_matrix.rowConst(face_list[l]);
- if (face_cells.length>1) {
+ const auto& face_cells = face_to_cell_matrix[face_list[l]];
+ if (face_cells.size()>1) {
std::cerr << "internal faces cannot be used to define mesh boundaries\n";
std::exit(1);
}
@@ -49,15 +49,14 @@ class MeshNodeBoundary
// not enough but should reduce significantly the number of resizing
node_ids.reserve(dimension*face_list.size());
const auto& face_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::face,
- ItemType::node);
+ = mesh.connectivity().faceToNodeMatrix();
for (size_t l=0; l<face_list.size(); ++l) {
const size_t face_number = face_list[l];
- const auto& face_nodes = face_to_node_matrix.rowConst(face_number);
+ const auto& face_nodes = face_to_node_matrix[face_number];
- for (size_t r=0; r<face_nodes.length; ++r) {
- node_ids.push_back(face_nodes(r));
+ for (size_t r=0; r<face_nodes.size(); ++r) {
+ node_ids.push_back(face_nodes[r]);
}
}
std::sort(node_ids.begin(), node_ids.end());
@@ -321,19 +320,17 @@ _getOutgoingNormal(const MeshType& mesh)
const NodeValue<const R>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[NodeId{j0_nodes(r)}]-xr[r0], normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -358,19 +355,17 @@ _getOutgoingNormal(const MeshType& mesh)
const NodeValue<const R2>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[NodeId{j0_nodes(r)}]-xr[r0], normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
@@ -395,19 +390,17 @@ _getOutgoingNormal(const MeshType& mesh)
const NodeValue<const R3>& xr = mesh.xr();
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
const auto& node_to_cell_matrix
- = mesh.connectivity().getMatrix(ItemType::node,
- ItemType::cell);
+ = mesh.connectivity().nodeToCellMatrix();
const NodeId r0 = m_node_list[0];
- const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
+ const CellId j0 = node_to_cell_matrix[r0][0];
+ const auto& j0_nodes = cell_to_node_matrix[j0];
double max_height = 0;
- for (size_t r=0; r<j0_nodes.length; ++r) {
- const double height = (xr[NodeId{j0_nodes(r)}]-xr[r0], normal);
+ for (size_t r=0; r<j0_nodes.size(); ++r) {
+ const double height = (xr[j0_nodes[r]]-xr[r0], normal);
if (std::abs(height) > std::abs(max_height)) {
max_height = height;
}
diff --git a/src/mesh/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
index 1c9a46ac0..3853038a1 100644
--- a/src/mesh/SubItemValuePerItem.hpp
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -213,7 +213,7 @@ class SubItemValuePerItem<DataType,
"Cannot allocate SubItemValuePerItem of const data: only view is supported"); ;
ConnectivityMatrix connectivity_matrix
- = connectivity.getMatrix(item_type, sub_item_type);
+ = connectivity._getMatrix(item_type, sub_item_type);
m_host_row_map = connectivity_matrix.rowsMap();
m_values = Array<std::remove_const_t<DataType>>(connectivity_matrix.numEntries());
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index f39fd3415..276bad33d 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -73,13 +73,12 @@ class VTKWriter
fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
const auto& cell_to_node_matrix
- = mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = mesh.connectivity().cellToNodeMatrix();
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- for (unsigned short r=0; r<cell_nodes.length; ++r) {
- fout << cell_nodes(r) << ' ';
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ for (unsigned short r=0; r<cell_nodes.size(); ++r) {
+ fout << cell_nodes[r] << ' ';
}
}
fout << '\n';
@@ -88,9 +87,9 @@ class VTKWriter
fout << "<DataArray type=\"UInt32\" Name=\"offsets\" NumberOfComponents=\"1\" format=\"ascii\">\n";
{
unsigned int offset=0;
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- offset += cell_nodes.length;
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ offset += cell_nodes.size();
fout << offset << ' ';
}
}
@@ -98,9 +97,9 @@ class VTKWriter
fout << "</DataArray>\n";
fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
- for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- switch (cell_nodes.length) {
+ for (CellId j=0; j<mesh.numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ switch (cell_nodes.size()) {
case 2: {
fout << "3 ";
break;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index e6d8d289d..64c3967cc 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -64,19 +64,18 @@ class AcousticSolver
const NodeValue<const Rdd>
computeAr(const NodeValuePerCell<const Rdd>& Ajr) {
const auto& node_to_cell_matrix
- = m_connectivity.getMatrix(ItemType::node,
- ItemType::cell);
+ = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
Rdd sum = zero;
- const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
+ const auto& node_to_cell = node_to_cell_matrix[r];//.rowConst(r);
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
- for (size_t j=0; j<node_to_cell.length; ++j) {
- const CellId J = node_to_cell(j);
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
sum += Ajr(J,R);
}
@@ -93,19 +92,18 @@ class AcousticSolver
const CellValue<const Rd>& uj,
const CellValue<const double>& pj) {
const auto& node_to_cell_matrix
- = m_connectivity.getMatrix(ItemType::node,
- ItemType::cell);
+ = m_connectivity.nodeToCellMatrix();
const auto& node_local_numbers_in_their_cells
= m_connectivity.nodeLocalNumbersInTheirCells();
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const NodeId& r) {
Rd& br = m_br[r];
br = zero;
- const auto& node_to_cell = node_to_cell_matrix.rowConst(r);
+ const auto& node_to_cell = node_to_cell_matrix[r];
const auto& node_local_number_in_its_cells
= node_local_numbers_in_their_cells.itemValues(r);
- for (size_t j=0; j<node_to_cell.length; ++j) {
- const CellId J = node_to_cell(j);
+ for (size_t j=0; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
const unsigned int R = node_local_number_in_its_cells[j];
br += Ajr(J,R)*uj[J] + pj[J]*Cjr(J,R);
}
@@ -179,14 +177,13 @@ class AcousticSolver
const CellValue<const double>& pj)
{
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
- for (size_t r=0; r<cell_nodes.length; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[NodeId{cell_nodes(r)}])+pj[j]*Cjr(j,r);
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
+ m_Fjr(j,r) = Ajr(j,r)*(uj[j]-ur[cell_nodes[r]])+pj[j]*Cjr(j,r);
}
});
}
@@ -259,14 +256,13 @@ class AcousticSolver
{
const NodeValuePerCell<const double>& ljr = m_mesh_data.ljr();
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j){
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
double S = 0;
- for (size_t r=0; r<cell_nodes.length; ++r) {
+ for (size_t r=0; r<cell_nodes.size(); ++r) {
S += ljr(j,r);
}
m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
@@ -298,17 +294,16 @@ class AcousticSolver
const NodeValuePerCell<Rd>& Fjr = m_Fjr;
const NodeValue<const Rd> ur = m_ur;
const auto& cell_to_node_matrix
- = m_mesh.connectivity().getMatrix(ItemType::cell,
- ItemType::node);
+ = m_mesh.connectivity().cellToNodeMatrix();
const CellValue<const double> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const CellId& j) {
- const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix[j];
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
- for (size_t R=0; R<cell_nodes.length; ++R) {
- const NodeId r=cell_nodes(R);
+ for (size_t R=0; R<cell_nodes.size(); ++R) {
+ const NodeId r=cell_nodes[R];
momentum_fluxes += Fjr(j,R);
energy_fluxes += (Fjr(j,R), ur[r]);
}
--
GitLab