From 5b7dfe1bf9bcb4186b854ba20782f15475f682c2 Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 1 Aug 2018 19:00:53 +0200
Subject: [PATCH] Begin encapsulation of item_to_item_matrix
- add access function to m_cell_to_node_matrix
- prepare array of ConnectivityMatrix and item type to item id mechanism
---
src/mesh/Connectivity.hpp | 17 +++++++--
src/mesh/MeshData.hpp | 16 ++++++---
src/mesh/MeshNodeBoundary.hpp | 15 ++++----
src/mesh/TypeOfItem.hpp | 65 +++++++++++++++++++++++++++++++++++
src/output/VTKWriter.hpp | 7 ++--
src/scheme/AcousticSolver.hpp | 9 +++--
6 files changed, 109 insertions(+), 20 deletions(-)
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 62710a761..d0448615b 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -224,12 +224,20 @@ template <size_t Dimension>
class Connectivity final
: public IConnectivity
{
-public:
+ private:
+ constexpr static auto& itemId = ItemId<Dimension>::itemId;
+
+ public:
static constexpr size_t dimension = Dimension;
- static constexpr bool has_edges = (dimension>2);
- static constexpr bool has_face = (dimension>1);
+ private:
ConnectivityMatrix m_cell_to_node_matrix;
+ public:
+ ConnectivityMatrix cellToNodeMatrix() const
+ {
+ return m_cell_to_node_matrix;
+ }
+
NodeValuePerCell<unsigned short> m_cell_to_node_local_cell;
ConnectivityMatrix m_cell_to_face_matrix;
@@ -250,6 +258,9 @@ public:
const TypeOfItem& item_type_1) const final;
private:
+
+ ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
+
ConnectivityComputer m_connectivity_computer;
std::vector<RefFaceList> m_ref_face_list;
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 49513c054..a3b1fb180 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -35,17 +35,20 @@ class MeshData
if constexpr (dimension == 1) {
const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
m_xj[j] = 0.5*(xr[cell_nodes(0)]+xr[cell_nodes(1)]);
});
} else {
const Kokkos::View<const Rd*> xr = m_mesh.xr();
const Kokkos::View<const double*>& inv_cell_nb_nodes
= m_mesh.connectivity().invCellNbNodes();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
Rd X = zero;
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (size_t R=0; R<cell_nodes.length; ++R) {
X += xr[cell_nodes(R)];
}
@@ -58,10 +61,11 @@ class MeshData
void _updateVolume()
{
const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
double sum_cjr_xr = 0;
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (size_t R=0; R<cell_nodes.length; ++R) {
sum_cjr_xr += (xr[cell_nodes(R)], m_Cjr(j,R));
@@ -77,9 +81,10 @@ class MeshData
}
else if constexpr (dimension == 2) {
const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(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;
@@ -128,8 +133,9 @@ class MeshData
m_Cjr[jr] = zero;
});
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
const auto& cell_faces = m_mesh.connectivity().m_cell_to_face_matrix.rowConst(j);
const auto& cell_face_is_reversed = m_mesh.connectivity().m_cell_face_is_reversed.itemValues(j);
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index b1cc055bd..a2b70d07c 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -306,15 +306,16 @@ _getOutgoingNormal(const MeshType& mesh)
const R normal = this->_getNormal(mesh);
const Kokkos::View<const R*>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const size_t r0 = m_node_list[0];
const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
const double height = (xr(j0_nodes(r))-xr(r0), normal);
if (std::abs(height) > std::abs(max_height)) {
- max_height = height;
+ max_height = height;
}
}
if (max_height > 0) {
@@ -336,15 +337,16 @@ _getOutgoingNormal(const MeshType& mesh)
const R2 normal = this->_getNormal(mesh);
const Kokkos::View<const R2*>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const size_t r0 = m_node_list[0];
const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
const double height = (xr(j0_nodes(r))-xr(r0), normal);
if (std::abs(height) > std::abs(max_height)) {
- max_height = height;
+ max_height = height;
}
}
if (max_height > 0) {
@@ -366,15 +368,16 @@ _getOutgoingNormal(const MeshType& mesh)
const R3 normal = this->_getNormal(mesh);
const Kokkos::View<const R3*>& xr = mesh.xr();
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
const size_t r0 = m_node_list[0];
const size_t j0 = mesh.connectivity().m_node_to_cell_matrix.rowConst(r0)(0);
- const auto& j0_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
for (size_t r=0; r<j0_nodes.length; ++r) {
const double height = (xr(j0_nodes(r))-xr(r0), normal);
if (std::abs(height) > std::abs(max_height)) {
- max_height = height;
+ max_height = height;
}
}
if (max_height > 0) {
diff --git a/src/mesh/TypeOfItem.hpp b/src/mesh/TypeOfItem.hpp
index 4b1b7d6a0..ccf4ef4f8 100644
--- a/src/mesh/TypeOfItem.hpp
+++ b/src/mesh/TypeOfItem.hpp
@@ -1,6 +1,8 @@
#ifndef TYPE_OF_ITEM_HPP
#define TYPE_OF_ITEM_HPP
+#include <utility>
+
enum class TypeOfItem {
node = 0,
edge = 1,
@@ -8,4 +10,67 @@ enum class TypeOfItem {
cell = 3
};
+template <size_t Dimension>
+struct ItemId {};
+
+template <>
+struct ItemId<1>
+{
+ static constexpr size_t itemId(const TypeOfItem& item_type) {
+ switch(item_type) {
+ case TypeOfItem::cell: {
+ return 0;
+ }
+ case TypeOfItem::face:
+ case TypeOfItem::edge:
+ case TypeOfItem::node: {
+ // in 1d, faces, edges and nodes are the same
+ return 1;
+ }
+ }
+ }
+};
+
+template <>
+struct ItemId<2>
+{
+ static constexpr size_t itemId(const TypeOfItem& item_type) {
+ switch(item_type) {
+ case TypeOfItem::cell: {
+ return 0;
+ }
+ case TypeOfItem::face:
+ case TypeOfItem::edge: {
+ // in 2d, faces and edges are the same
+ return 1;
+ }
+ case TypeOfItem::node: {
+ return 2;
+ }
+ }
+ }
+};
+
+template <>
+struct ItemId<3>
+{
+ static constexpr size_t itemId(const TypeOfItem& item_type) {
+ switch(item_type) {
+ case TypeOfItem::cell: {
+ return 0;
+ }
+ case TypeOfItem::face: {
+ return 1;
+ }
+ case TypeOfItem::edge: {
+ // in 2d, faces and edges are the same
+ return 2;
+ }
+ case TypeOfItem::node: {
+ return 3;
+ }
+ }
+ }
+};
+
#endif // TYPE_OF_ITEM_HPP
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index f297325bd..d4ea1f33a 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -69,8 +69,9 @@ class VTKWriter
fout << "<Cells>\n";
fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+ const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(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) << ' ';
}
@@ -82,7 +83,7 @@ class VTKWriter
{
unsigned int offset=0;
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
offset += cell_nodes.length;
fout << offset << ' ';
}
@@ -92,7 +93,7 @@ class VTKWriter
fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- const auto& cell_nodes = mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
switch (cell_nodes.length) {
case 2: {
fout << "3 ";
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 416d6c642..2ad5c7203 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -196,8 +196,9 @@ class AcousticSolver
const Kokkos::View<const Rd*>& uj,
const Kokkos::View<const double*>& pj)
{
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
for (size_t r=0; r<cell_nodes.length; ++r) {
m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(r)))+pj(j)*Cjr(j,r);
@@ -271,9 +272,10 @@ class AcousticSolver
const Kokkos::View<const double*>& cj) const
{
const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
double S = 0;
for (size_t r=0; r<cell_nodes.length; ++r) {
@@ -312,10 +314,11 @@ class AcousticSolver
const NodeValuePerCell<Rd>& Fjr = m_Fjr;
const Kokkos::View<const Rd*> ur = m_ur;
+ const auto& cell_to_node_matrix = m_mesh.connectivity().cellToNodeMatrix();
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- const auto& cell_nodes = m_mesh.connectivity().m_cell_to_node_matrix.rowConst(j);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
Rd momentum_fluxes = zero;
double energy_fluxes = 0;
--
GitLab