diff --git a/src/main.cpp b/src/main.cpp
index f277ce412860ccff5b88c256a2b552a8443e3e73..b1a3decd13b36a4b738c49c9815a0c4bfb327d5f 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -11,9 +11,7 @@
// #include <AcousticSolverClass.hpp>
// #include <AcousticSolverTest.hpp>
-#include <Connectivity1D.hpp>
-#include <Connectivity2D.hpp>
-#include <Connectivity3D.hpp>
+#include <Connectivity.hpp>
#include <Mesh.hpp>
#include <BoundaryCondition.hpp>
@@ -197,7 +195,6 @@ int main(int argc, char *argv[])
typedef TinyVector<MeshType::dimension> Rd;
const Kokkos::View<const double*> Vj = mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.2;
double t=0;
@@ -304,10 +301,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- typedef TinyVector<MeshType::dimension> Rd;
-
const Kokkos::View<const double*> Vj = mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.2;
double t=0;
@@ -403,10 +397,7 @@ int main(int argc, char *argv[])
AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
- typedef TinyVector<MeshType::dimension> Rd;
-
const Kokkos::View<const double*> Vj = mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
const double tmax=0.2;
double t=0;
@@ -431,7 +422,6 @@ int main(int argc, char *argv[])
dt=tmax-t;
}
acoustic_solver.computeNextStep(t,dt, unknowns);
-
block_eos.updatePandCFromRhoE();
t += dt;
@@ -450,90 +440,8 @@ int main(int argc, char *argv[])
std::cout << "* " << rang::fgB::red << "Could not be uglier!" << rang::fg::reset << " (" << __FILE__ << ':' << __LINE__ << ")\n";
} else {
- // class for acoustic solver test
- Kokkos::Timer timer;
- timer.reset();
- std::shared_ptr<Connectivity1D>connectivity( new Connectivity1D(number));
- typedef Mesh<Connectivity1D> MeshType;
- typedef MeshData<MeshType> MeshDataType;
- typedef FiniteVolumesEulerUnknowns<MeshDataType> UnknownsType;
-
- MeshType mesh(connectivity);
- MeshDataType mesh_data(mesh);
-
- std::vector<BoundaryConditionHandler> bc_list;
- { // quite dirty!
- // SymmetryBoundaryCondition<MeshType::dimension>* sym_bc0
- // = new SymmetryBoundaryCondition<MeshType::dimension>(std::vector<unsigned int>({0u}),
- // TinyVector<1>(-1));
- // std::shared_ptr<SymmetryBoundaryCondition<1>> bc0(sym_bc0);
- // bc_list.push_back(BoundaryConditionHandler(bc0));
- PressureBoundaryCondition* pres_bc0
- = new PressureBoundaryCondition(1,
- std::vector<unsigned int>({static_cast<unsigned int>(0u)}));
- std::shared_ptr<PressureBoundaryCondition> bc0(pres_bc0);
- bc_list.push_back(BoundaryConditionHandler(bc0));
-
- PressureBoundaryCondition* pres_bc1
- = new PressureBoundaryCondition(0.1,
- std::vector<unsigned int>({static_cast<unsigned int>(mesh.numberOfCells())}));
- std::shared_ptr<PressureBoundaryCondition> bc1(pres_bc1);
- bc_list.push_back(BoundaryConditionHandler(bc1));
- }
-
- UnknownsType unknowns(mesh_data);
-
- unknowns.initializeSod();
-
- AcousticSolver<MeshDataType> acoustic_solver(mesh_data, unknowns, bc_list);
-
- typedef TinyVector<MeshType::dimension> Rd;
-
- const Kokkos::View<const double*> Vj = mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = mesh_data.Cjr();
-
- const double tmax=0.2;
- double t=0;
-
- int itermax=std::numeric_limits<int>::max();
- int iteration=0;
-
- Kokkos::View<double*> rhoj = unknowns.rhoj();
- Kokkos::View<double*> ej = unknowns.ej();
- Kokkos::View<double*> pj = unknowns.pj();
- Kokkos::View<double*> gammaj = unknowns.gammaj();
- Kokkos::View<double*> cj = unknowns.cj();
-
- BlockPerfectGas block_eos(rhoj, ej, pj, gammaj, cj);
-
- while((t<tmax) and (iteration<itermax)) {
- double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
- if (t+dt>tmax) {
- dt=tmax-t;
- }
-
- acoustic_solver.computeNextStep(t,dt, unknowns);
-
- block_eos.updatePandCFromRhoE();
-
- t += dt;
- ++iteration;
- }
-
- std::cout << "* " << rang::style::underline << "Final time" << rang::style::reset
- << ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
-
- method_cost_map["AcousticSolverWithMesh"] = timer.seconds();
-
- { // gnuplot output for density
- const Kokkos::View<const Rd*> xj = mesh_data.xj();
- const Kokkos::View<const double*> rhoj = unknowns.rhoj();
- std::ofstream fout("rho");
- for (size_t j=0; j<mesh.numberOfCells(); ++j) {
- fout << xj[j][0] << ' ' << rhoj[j] << '\n';
- }
- }
-
+ std::cerr << "Connectivity1D defined by number of nodes no more implemented\n";
+ std::exit(0);
}
}
catch (const AssertError& error) {
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index f07bc30ba6e13173157a7ddc0b0664ea0920bf60..c28b0643fd8bbf1b2677512423d2185d9a7ff970 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -5,6 +5,8 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR})
add_library(
PastisMesh
+ Connectivity.cpp
+ ConnectivityComputer.cpp
GmshReader.cpp)
#include_directories(${PASTIS_SOURCE_DIR}/utils)
diff --git a/src/mesh/CellType.hpp b/src/mesh/CellType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..1f444323b5f5073b9aa09d020cc1fc962b1182a5
--- /dev/null
+++ b/src/mesh/CellType.hpp
@@ -0,0 +1,17 @@
+#ifndef CELL_TYPE_HPP
+#define CELL_TYPE_HPP
+
+enum class CellType : unsigned short
+{
+ Line,
+
+ Triangle,
+ Quadrangle,
+
+ Tetrahedron,
+ Pyramid,
+ Prism,
+ Hexahedron
+};
+
+#endif // CELL_TYPE_HPP
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..eb5fb87b36947a5e4778f2e38161e6f35cb3f92d
--- /dev/null
+++ b/src/mesh/Connectivity.cpp
@@ -0,0 +1,261 @@
+#include <Connectivity.hpp>
+#include <map>
+
+template<>
+void Connectivity<3>::_computeCellFaceAndFaceNodeConnectivities()
+{
+ typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
+
+ const auto& cell_to_node_matrix
+ = this->getMatrix(ItemType::cell, ItemType::node);
+
+ Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", this->numberOfCells());
+ std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
+ for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+
+ switch (m_cell_type[j]) {
+ case CellType::Tetrahedron: {
+ cell_nb_faces[j] = 4;
+ // face 0
+ Face f0({cell_nodes(1),
+ cell_nodes(2),
+ cell_nodes(3)});
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes(0),
+ cell_nodes(3),
+ cell_nodes(2)});
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes(0),
+ cell_nodes(1),
+ cell_nodes(3)});
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+ // face 3
+ Face f3({cell_nodes(0),
+ cell_nodes(2),
+ cell_nodes(1)});
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ break;
+ }
+ case CellType::Hexahedron: {
+ // face 0
+ Face f0({cell_nodes(3),
+ cell_nodes(2),
+ cell_nodes(1),
+ cell_nodes(0)});
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+
+ // face 1
+ Face f1({cell_nodes(4),
+ cell_nodes(5),
+ cell_nodes(6),
+ cell_nodes(7)});
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+
+ // face 2
+ Face f2({cell_nodes(0),
+ cell_nodes(4),
+ cell_nodes(7),
+ cell_nodes(3)});
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+
+ // face 3
+ Face f3({cell_nodes(1),
+ cell_nodes(2),
+ cell_nodes(6),
+ cell_nodes(5)});
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+
+ // face 4
+ Face f4({cell_nodes(0),
+ cell_nodes(1),
+ cell_nodes(5),
+ cell_nodes(4)});
+ face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+
+ // face 5
+ Face f5({cell_nodes(3),
+ cell_nodes(7),
+ cell_nodes(6),
+ cell_nodes(2)});
+ face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+
+ cell_nb_faces[j] = 6;
+ break;
+ }
+ default: {
+ std::cerr << "unexpected cell type!\n";
+ std::exit(0);
+ }
+ }
+ }
+
+ {
+ auto& cell_to_face_matrix
+ = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::face)];
+ std::vector<std::vector<unsigned int>> cell_to_face_vector(this->numberOfCells());
+ for (size_t j=0; j<cell_to_face_vector.size(); ++j) {
+ cell_to_face_vector[j].resize(cell_nb_faces[j]);
+ }
+ int l=0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& cells_vector = face_cells_vector.second;
+ for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+ const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+ cell_to_face_vector[cell_number][cell_local_face] = l;
+ }
+ ++l;
+ }
+ cell_to_face_matrix = cell_to_face_vector;
+ }
+
+ FaceValuePerCell<bool> cell_face_is_reversed(*this);
+ {
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& cells_vector = face_cells_vector.second;
+ for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
+ const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
+ cell_face_is_reversed(cell_number, cell_local_face) = reversed;
+ }
+ }
+
+ m_cell_face_is_reversed = cell_face_is_reversed;
+ }
+
+ {
+ auto& face_to_node_matrix
+ = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::node)];
+
+ std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
+ int l=0;
+ for (const auto& face_info : face_cells_map) {
+ const Face& face = face_info.first;
+ face_to_node_vector[l] = face.nodeIdList();
+ ++l;
+ }
+ face_to_node_matrix = face_to_node_vector;
+ }
+
+ {
+ int l=0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const Face& face = face_cells_vector.first;
+ m_face_number_map[face] = l;
+ ++l;
+ }
+ }
+
+#warning check that the number of cell per faces is <=2
+}
+
+template<>
+void Connectivity<2>::_computeCellFaceAndFaceNodeConnectivities()
+{
+ const auto& cell_to_node_matrix
+ = this->getMatrix(ItemType::cell, ItemType::node);
+
+ // In 2D faces are simply define
+ typedef std::pair<unsigned int, unsigned short> CellFaceId;
+ std::map<Face, std::vector<CellFaceId>> face_cells_map;
+ for (unsigned int j=0; j<this->numberOfCells(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ for (unsigned short r=0; r<cell_nodes.length; ++r) {
+ unsigned int node0_id = cell_nodes(r);
+ unsigned int node1_id = cell_nodes((r+1)%cell_nodes.length);
+ if (node1_id<node0_id) {
+ std::swap(node0_id, node1_id);
+ }
+ face_cells_map[Face({node0_id, node1_id})].push_back(std::make_pair(j, r));
+ }
+ }
+
+ {
+ int l=0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const Face& face = face_cells_vector.first;
+ m_face_number_map[face] = l;
+ ++l;
+ }
+ }
+
+ {
+ std::vector<std::vector<unsigned int>> face_to_node_vector(face_cells_map.size());
+ int l=0;
+ for (const auto& face_info : face_cells_map) {
+ const Face& face = face_info.first;
+ face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
+ ++l;
+ }
+ auto& face_to_node_matrix
+ = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::node)];
+ face_to_node_matrix = face_to_node_vector;
+ }
+
+ {
+ std::vector<std::vector<unsigned int>> face_to_cell_vector(face_cells_map.size());
+ int l=0;
+ for (const auto& face_cells_vector : face_cells_map) {
+ const auto& [face, cell_info_vector] = face_cells_vector;
+ for (const auto& cell_info : cell_info_vector) {
+ face_to_cell_vector[l].push_back(cell_info.second);
+ }
+ ++l;
+ }
+ auto& face_to_cell_matrix
+ = m_item_to_item_matrix[itemId(ItemType::face)][itemId(ItemType::cell)];
+ face_to_cell_matrix = face_to_cell_vector;
+ }
+}
+
+
+template<size_t Dimension>
+Connectivity<Dimension>::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+ const std::vector<CellType>& cell_type_vector)
+{
+ Assert(cell_by_node_vector.size() == cell_type_vector.size());
+
+ auto& cell_to_node_matrix
+ = m_item_to_item_matrix[itemId(ItemType::cell)][itemId(ItemType::node)];
+ cell_to_node_matrix = cell_by_node_vector;
+
+ Assert(this->numberOfCells()>0);
+
+ {
+ Kokkos::View<CellType*> cell_type("cell_type", this->numberOfCells());
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ cell_type[j] = cell_type_vector[j];
+ });
+ m_cell_type = cell_type;
+ }
+ {
+ Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", this->numberOfCells());
+ Kokkos::parallel_for(this->numberOfCells(), KOKKOS_LAMBDA(const int& j){
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ inv_cell_nb_nodes[j] = 1./cell_nodes.length;
+ });
+ m_inv_cell_nb_nodes = inv_cell_nb_nodes;
+ }
+
+ if constexpr (Dimension>1) {
+ this->_computeCellFaceAndFaceNodeConnectivities();
+ }
+}
+
+
+template Connectivity1D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+ const std::vector<CellType>& cell_type_vector);
+
+template Connectivity2D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+ const std::vector<CellType>& cell_type_vector);
+
+template Connectivity3D::
+Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+ const std::vector<CellType>& cell_type_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0994fc9c84e4105c3f9de66b55231f3fd61ff8fe
--- /dev/null
+++ b/src/mesh/Connectivity.hpp
@@ -0,0 +1,498 @@
+#ifndef CONNECTIVITY_HPP
+#define CONNECTIVITY_HPP
+
+#include <PastisAssert.hpp>
+#include <TinyVector.hpp>
+
+#include <Kokkos_Core.hpp>
+
+#include <IConnectivity.hpp>
+
+#include <ConnectivityMatrix.hpp>
+#include <SubItemValuePerItem.hpp>
+#include <ConnectivityComputer.hpp>
+
+#include <vector>
+#include <unordered_map>
+#include <algorithm>
+
+#include <CellType.hpp>
+
+#include <RefId.hpp>
+#include <ItemType.hpp>
+#include <RefNodeList.hpp>
+#include <RefFaceList.hpp>
+
+#include <tuple>
+#include <algorithm>
+
+template <size_t Dimension>
+class Connectivity;
+
+template <size_t Dimension>
+class ConnectivityFace;
+
+template<>
+class ConnectivityFace<1>
+{
+ public:
+ friend struct Hash;
+ struct Hash
+ {
+ size_t operator()(const ConnectivityFace& f) const;
+ };
+};
+
+template<>
+class ConnectivityFace<2>
+{
+ public:
+ friend struct Hash;
+ struct Hash
+ {
+ size_t operator()(const ConnectivityFace& f) const {
+ size_t hash = 0;
+ hash ^= std::hash<unsigned int>()(f.m_node0_id);
+ hash ^= std::hash<unsigned int>()(f.m_node1_id) >> 1;
+ return hash;
+ }
+ };
+
+ unsigned int m_node0_id;
+ unsigned int m_node1_id;
+
+ friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
+ {
+ os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
+ return os;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ bool operator==(const ConnectivityFace& f) const
+ {
+ return ((m_node0_id == f.m_node0_id) and
+ (m_node1_id == f.m_node1_id));
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ bool operator<(const ConnectivityFace& f) const
+ {
+ return ((m_node0_id<f.m_node0_id) or
+ ((m_node0_id == f.m_node0_id) and
+ (m_node1_id<f.m_node1_id)));
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace& operator=(const ConnectivityFace&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace& operator=(ConnectivityFace&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
+ {
+ Assert(given_node_id_list.size()==2);
+#warning rework this dirty constructor
+ const auto& [min, max] = std::minmax(given_node_id_list[0], given_node_id_list[1]);
+ m_node0_id = min;
+ m_node1_id = max;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(const ConnectivityFace&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(ConnectivityFace&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ~ConnectivityFace() = default;
+};
+
+template <>
+class ConnectivityFace<3>
+{
+ private:
+ friend class Connectivity<3>;
+ friend struct Hash;
+ struct Hash
+ {
+ size_t operator()(const ConnectivityFace& f) const {
+ size_t hash = 0;
+ for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
+ hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
+ }
+ return hash;
+ }
+ };
+
+ bool m_reversed;
+ std::vector<unsigned int> m_node_id_list;
+
+ friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
+ {
+ for (auto id : f.m_node_id_list) {
+ std::cout << id << ' ';
+ }
+ return os;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const bool& reversed() const
+ {
+ return m_reversed;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const std::vector<unsigned int>& nodeIdList() const
+ {
+ return m_node_id_list;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
+ {
+ const auto min_id = std::min_element(node_list.begin(), node_list.end());
+ const int shift = std::distance(node_list.begin(), min_id);
+
+ std::vector<unsigned int> rotated_node_list(node_list.size());
+ if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
+ for (size_t i=0; i<node_list.size(); ++i) {
+ rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
+ m_reversed = true;
+ }
+ } else {
+ for (size_t i=0; i<node_list.size(); ++i) {
+ rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
+ }
+ }
+
+ return rotated_node_list;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
+ : m_reversed(false),
+ m_node_id_list(_sort(given_node_id_list))
+ {
+ ;
+ }
+
+ public:
+ bool operator==(const ConnectivityFace& f) const
+ {
+ if (m_node_id_list.size() == f.nodeIdList().size()) {
+ for (size_t j=0; j<m_node_id_list.size(); ++j) {
+ if (m_node_id_list[j] != f.nodeIdList()[j]) {
+ return false;
+ }
+ }
+ return true;
+ }
+ return false;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ bool operator<(const ConnectivityFace& f) const
+ {
+ const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
+ for (size_t i=0; i<min_nb_nodes; ++i) {
+ if (m_node_id_list[i] < f.m_node_id_list[i]) return true;
+ if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
+ }
+ return m_node_id_list.size() < f.m_node_id_list.size();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace& operator=(const ConnectivityFace&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace& operator=(ConnectivityFace&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(const ConnectivityFace&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace(ConnectivityFace&&) = default;
+
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityFace() = delete;
+
+ KOKKOS_INLINE_FUNCTION
+ ~ConnectivityFace() = default;
+};
+
+template <size_t Dimension>
+class Connectivity final
+ : public IConnectivity
+{
+ private:
+ constexpr static auto& itemId = ItemId<Dimension>::itemId;
+
+ public:
+ static constexpr size_t dimension = Dimension;
+
+ private:
+ ConnectivityMatrix m_item_to_item_matrix[Dimension+1][Dimension+1];
+ Kokkos::View<const CellType*> m_cell_type;
+
+ FaceValuePerCell<const bool> m_cell_face_is_reversed;
+
+ NodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
+ EdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
+ FaceValuePerCell<const unsigned short> m_cell_local_numbers_in_their_faces;
+
+ CellValuePerFace<const unsigned short> m_face_local_numbers_in_their_cells;
+ EdgeValuePerFace<const unsigned short> m_face_local_numbers_in_their_edges;
+ NodeValuePerFace<const unsigned short> m_face_local_numbers_in_their_nodes;
+
+ CellValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_cells;
+ FaceValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_faces;
+ NodeValuePerEdge<const unsigned short> m_edge_local_numbers_in_their_nodes;
+
+ CellValuePerNode<const unsigned short> m_node_local_numbers_in_their_cells;
+ EdgeValuePerNode<const unsigned short> m_node_local_numbers_in_their_edges;
+ FaceValuePerNode<const unsigned short> m_node_local_numbers_in_their_faces;
+
+ ConnectivityComputer m_connectivity_computer;
+
+ std::vector<RefFaceList> m_ref_face_list;
+ std::vector<RefNodeList> m_ref_node_list;
+
+ Kokkos::View<double*> m_inv_cell_nb_nodes;
+
+ using Face = ConnectivityFace<Dimension>;
+
+ std::unordered_map<Face, unsigned int, typename Face::Hash> m_face_number_map;
+
+ void _computeCellFaceAndFaceNodeConnectivities();
+
+ template <typename SubItemValuePerItemType>
+ KOKKOS_INLINE_FUNCTION
+ const SubItemValuePerItemType&
+ _lazzyBuildItemNumberInTheirChild(const SubItemValuePerItemType& sub_item_value_per_item) const
+ {
+ if (not sub_item_value_per_item.isBuilt()) {
+ const_cast<SubItemValuePerItemType&>(sub_item_value_per_item)
+ = m_connectivity_computer
+ . computeLocalItemNumberInChildItem<SubItemValuePerItemType::item_t,
+ SubItemValuePerItemType::sub_item_t>(*this);
+ }
+ return sub_item_value_per_item;
+ }
+
+ public:
+
+ KOKKOS_INLINE_FUNCTION
+ const bool& isConnectivityMatrixBuilt(const ItemType& item_type_0,
+ const ItemType& item_type_1) const
+ {
+ const ConnectivityMatrix& connectivity_matrix
+ = m_item_to_item_matrix[itemId(item_type_0)][itemId(item_type_1)];
+ return connectivity_matrix.isBuilt();
+ }
+
+ 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[itemId(item_type_0)][itemId(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;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& cellFaceIsReversed() const
+ {
+ static_assert(dimension>1, "reversed faces makes no sense in dimension 1");
+ return m_cell_face_is_reversed;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& cellLocalNumbersInTheirNodes() const
+ {
+ return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_nodes);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& cellLocalNumbersInTheirEdges() const
+ {
+ if constexpr (dimension>2) {
+ return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_edges);
+ } else {
+ return cellLocalNumbersInTheirFaces();
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& cellLocalNumbersInTheirFaces() const
+ {
+ if constexpr (dimension>1) {
+ return _lazzyBuildItemNumberInTheirChild(m_cell_local_numbers_in_their_faces);
+ } else {
+ return cellLocalNumbersInTheirNodes();
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& faceLocalNumbersInTheirCells() const
+ {
+ if constexpr(dimension>1) {
+ return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_cells);
+ } else {
+ return nodeLocalNumbersInTheirCells();
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& faceLocalNumbersInTheirEdges() const
+ {
+ static_assert(dimension>2,"this function has no meaning in 1d or 2d");
+ return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_edges);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& faceLocalNumbersInTheirNodes() const
+ {
+ static_assert(dimension>1,"this function has no meaning in 1d");
+ return _lazzyBuildItemNumberInTheirChild(m_face_local_numbers_in_their_nodes);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& edgeLocalNumbersInTheirCells() const
+ {
+ if constexpr (dimension>2) {
+ return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_cells);
+ } else {
+ return faceLocalNumbersInTheirCells();
+ }
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& edgeLocalNumbersInTheirFaces() const
+ {
+ static_assert(dimension>2, "this function has no meaning in 1d or 2d");
+ return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_faces);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& edgeLocalNumbersInTheirNodes() const
+ {
+ static_assert(dimension>2, "this function has no meaning in 1d or 2d");
+ return _lazzyBuildItemNumberInTheirChild(m_edge_local_numbers_in_their_nodes);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& nodeLocalNumbersInTheirCells() const
+ {
+ return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_cells);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& nodeLocalNumbersInTheirEdges() const
+ {
+ static_assert(dimension>2, "this function has no meaning in 1d or 2d");
+ return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_edges);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto& nodeLocalNumbersInTheirFaces() const
+ {
+ static_assert(dimension>1,"this function has no meaning in 1d");
+ return _lazzyBuildItemNumberInTheirChild(m_node_local_numbers_in_their_faces);
+ }
+
+ void addRefFaceList(const RefFaceList& ref_face_list)
+ {
+ m_ref_face_list.push_back(ref_face_list);
+ }
+
+ size_t numberOfRefFaceList() const
+ {
+ return m_ref_face_list.size();
+ }
+
+ const RefFaceList& refFaceList(const size_t& i) const
+ {
+ return m_ref_face_list[i];
+ }
+
+ void addRefNodeList(const RefNodeList& ref_node_list)
+ {
+ m_ref_node_list.push_back(ref_node_list);
+ }
+
+ size_t numberOfRefNodeList() const
+ {
+ return m_ref_node_list.size();
+ }
+
+ const RefNodeList& refNodeList(const size_t& i) const
+ {
+ return m_ref_node_list[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfNodes() const
+ {
+ const auto& node_to_cell_matrix
+ = this->getMatrix(ItemType::node,ItemType::cell);
+ return node_to_cell_matrix.numRows();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfFaces() const
+ {
+ const auto& face_to_node_matrix
+ = this->getMatrix(ItemType::face,ItemType::cell);
+ return face_to_node_matrix.numRows();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfCells() const
+ {
+ const auto& cell_to_node_matrix
+ = this->getMatrix(ItemType::cell,ItemType::node);
+ return cell_to_node_matrix.numRows();
+ }
+
+ const Kokkos::View<const double*> invCellNbNodes() const
+ {
+#warning add calculation on demand when variables will be defined
+ return m_inv_cell_nb_nodes;
+ }
+
+ unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const
+ {
+ const Face face(face_nodes);
+ auto i_face = m_face_number_map.find(face);
+ if (i_face == m_face_number_map.end()) {
+ std::cerr << "Face " << face << " not found!\n";
+ throw std::exception();
+ std::exit(0);
+ }
+ return i_face->second;
+ }
+
+ Connectivity(const Connectivity&) = delete;
+
+ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
+ const std::vector<CellType>& cell_type_vector);
+
+ ~Connectivity()
+ {
+ ;
+ }
+};
+
+using Connectivity3D = Connectivity<3>;
+using Connectivity2D = Connectivity<2>;
+using Connectivity1D = Connectivity<1>;
+
+#endif // CONNECTIVITY_HPP
diff --git a/src/mesh/Connectivity1D.hpp b/src/mesh/Connectivity1D.hpp
deleted file mode 100644
index d0422f7b25a7b1df91e444f9bea0b9e8ab04a602..0000000000000000000000000000000000000000
--- a/src/mesh/Connectivity1D.hpp
+++ /dev/null
@@ -1,211 +0,0 @@
-#ifndef CONNECTIVITY_1D_HPP
-#define CONNECTIVITY_1D_HPP
-
-#include <Kokkos_Core.hpp>
-#include <PastisAssert.hpp>
-
-#include <TinyVector.hpp>
-#include <ConnectivityUtils.hpp>
-
-#include <RefId.hpp>
-#include <RefNodeList.hpp>
-
-class Connectivity1D
-{
-public:
- static constexpr size_t dimension = 1;
-
-private:
- std::vector<RefNodeList> m_ref_node_list;
-
- size_t m_number_of_nodes;
- const size_t& m_number_of_faces = m_number_of_nodes;
- const size_t m_number_of_cells;
-
- const Kokkos::View<const unsigned int**> m_cell_nodes;
- const Kokkos::View<const unsigned int**>& m_cell_faces = m_cell_nodes;
-
- const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
- Kokkos::View<double*> m_inv_cell_nb_nodes;
- const Kokkos::View<const unsigned short*>& m_cell_nb_faces = m_cell_nb_nodes;
-
- Kokkos::View<const unsigned short*> m_node_nb_cells;
- const Kokkos::View<const unsigned short*>& m_face_nb_cells = m_node_nb_cells;
-
- Kokkos::View<const unsigned int**> m_node_cells;
- const Kokkos::View<const unsigned int**>& m_face_cells = m_node_cells;
-
- Kokkos::View<const unsigned short**> m_node_cell_local_node;
- const Kokkos::View<const unsigned short**>& m_face_cell_local_face = m_node_cell_local_node;
-
- size_t m_max_nb_node_per_cell;
-
- const Kokkos::View<const unsigned int**>
- _buildCellNodes(const size_t& number_of_cells)
- {
- Kokkos::View<unsigned int*[2]> cell_nodes("cell_nodes", number_of_cells);
-
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) {
- cell_nodes(j,0) = j;
- cell_nodes(j,1) = j+1;
- });
-
- return cell_nodes;
- }
-
-
- const Kokkos::View<const unsigned short*>
- _buildCellNbNodes(const size_t& number_of_cells)
- {
- Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", number_of_cells);
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) {
- cell_nb_nodes[j] = 2;
- });
-
- return cell_nb_nodes;
- }
-
-public:
- void addRefNodeList(const RefNodeList& ref_node_list)
- {
- m_ref_node_list.push_back(ref_node_list);
- }
-
- size_t numberOfRefNodeList() const
- {
- return m_ref_node_list.size();
- }
-
- const RefNodeList& refNodeList(const size_t& i) const
- {
- return m_ref_node_list[i];
- }
-
- const size_t& numberOfNodes() const
- {
- return m_number_of_nodes;
- }
-
- const size_t& numberOfFaces() const
- {
- return m_number_of_faces;
- }
-
- const size_t& numberOfCells() const
- {
- return m_number_of_cells;
- }
-
- const size_t& maxNbNodePerCell() const
- {
- return m_max_nb_node_per_cell;
- }
-
- const Kokkos::View<const unsigned int**> cellNodes() const
- {
- return m_cell_nodes;
- }
-
- const Kokkos::View<const unsigned int**> cellFaces() const
- {
- return m_cell_faces;
- }
-
- const Kokkos::View<const unsigned short*> nodeNbCells() const
- {
- return m_node_nb_cells;
- }
-
- const Kokkos::View<const unsigned short*> cellNbNodes() const
- {
- return m_cell_nb_nodes;
- }
-
- const Kokkos::View<const double*> invCellNbNodes() const
- {
- return m_inv_cell_nb_nodes;
- }
-
- const Kokkos::View<const unsigned short*> cellNbFaces() const
- {
- return m_cell_nb_faces;
- }
-
- const Kokkos::View<const unsigned short*> faceNbCells() const
- {
- return m_face_nb_cells;
- }
-
- const Kokkos::View<const unsigned int**> nodeCells() const
- {
- return m_node_cells;
- }
-
- const Kokkos::View<const unsigned int**> faceCells() const
- {
- return m_face_cells;
- }
-
- const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
- {
- return m_node_cell_local_node;
- }
-
- const Kokkos::View<const unsigned short**> faceCellLocalFace() const
- {
- return m_face_cell_local_face;
- }
-
- Connectivity1D(const Connectivity1D&) = delete;
-
- Connectivity1D(const size_t& number_of_cells)
- : m_number_of_cells (number_of_cells),
- m_cell_nodes (_buildCellNodes(number_of_cells)),
- m_cell_nb_nodes (_buildCellNbNodes(number_of_cells)),
- m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells)
- {
- Assert(number_of_cells>0);
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) {
- m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
- });
-
- ConnectivityUtils utils;
- utils.computeNodeCellConnectivity(m_cell_nodes,
- m_cell_nb_nodes,
- m_number_of_cells,
- m_max_nb_node_per_cell,
- m_number_of_nodes,
- m_node_nb_cells,
- m_node_cells,
- m_node_cell_local_node);
- }
-
- Connectivity1D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
- const Kokkos::View<const unsigned int**> cell_nodes)
- : m_number_of_cells (cell_nb_nodes.extent(0)),
- m_cell_nodes (cell_nodes),
- m_cell_nb_nodes (cell_nb_nodes),
- m_inv_cell_nb_nodes ("inv_cell_nb_nodes", m_number_of_cells)
- {
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const size_t& j) {
- m_inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
- });
-
- ConnectivityUtils utils;
- utils.computeNodeCellConnectivity(m_cell_nodes,
- m_cell_nb_nodes,
- m_number_of_cells,
- m_max_nb_node_per_cell,
- m_number_of_nodes,
- m_node_nb_cells,
- m_node_cells,
- m_node_cell_local_node);
- }
-
- ~Connectivity1D()
- {
- ;
- }
-};
-
-#endif // CONNECTIVITY_1D_HPP
diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
deleted file mode 100644
index bf87b5b9ab2f69306ec4dbc68958927f7df6df6f..0000000000000000000000000000000000000000
--- a/src/mesh/Connectivity2D.hpp
+++ /dev/null
@@ -1,352 +0,0 @@
-#ifndef CONNECTIVITY_2D_HPP
-#define CONNECTIVITY_2D_HPP
-
-#include <Kokkos_Core.hpp>
-#include <PastisAssert.hpp>
-#include <TinyVector.hpp>
-
-#include <ConnectivityUtils.hpp>
-#include <vector>
-#include <map>
-#include <algorithm>
-
-#include <RefId.hpp>
-#include <RefNodeList.hpp>
-#include <RefFaceList.hpp>
-
-class Connectivity2D
-{
- public:
- static constexpr size_t dimension = 2;
-
- private:
- std::vector<RefFaceList> m_ref_face_list;
- std::vector<RefNodeList> m_ref_node_list;
-
- const size_t m_number_of_cells;
- size_t m_number_of_faces;
- size_t m_number_of_nodes;
-
- const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
- const Kokkos::View<const unsigned int**> m_cell_nodes;
- Kokkos::View<double*> m_inv_cell_nb_nodes;
-
- Kokkos::View<const unsigned short*> m_cell_nb_faces;
- Kokkos::View<unsigned int**> m_cell_faces;
-
- Kokkos::View<const unsigned short*> m_node_nb_cells;
- Kokkos::View<const unsigned int**> m_node_cells;
- Kokkos::View<const unsigned short**> m_node_cell_local_node;
-
- Kokkos::View<unsigned short*> m_face_nb_cells;
- Kokkos::View<unsigned int**> m_face_cells;
- Kokkos::View<unsigned short**> m_face_cell_local_face;
-
- Kokkos::View<unsigned short*> m_face_nb_nodes;
- Kokkos::View<unsigned int**> m_face_nodes;
- Kokkos::View<unsigned short**> m_face_node_local_face;
-
- size_t m_max_nb_node_per_cell;
-
- struct Face
- {
- const unsigned int m_node0_id;
- const unsigned int m_node1_id;
-
- KOKKOS_INLINE_FUNCTION
- bool operator<(const Face& f) const
- {
- return ((m_node0_id<f.m_node0_id) or
- ((m_node0_id == f.m_node0_id) and
- (m_node1_id<f.m_node1_id)));
- }
-
- KOKKOS_INLINE_FUNCTION
- Face& operator=(const Face&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face& operator=(Face&&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(const Face&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(Face&&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(unsigned int node0_id,
- unsigned int node1_id)
- : m_node0_id(node0_id),
- m_node1_id(node1_id)
- {
- ;
- }
-
- KOKKOS_INLINE_FUNCTION
- ~Face() = default;
- };
-
- void _computeFaceCellConnectivities()
- {
- // In 2D faces are simply define
- typedef std::pair<unsigned int, unsigned short> CellFaceId;
- std::map<Face, std::vector<CellFaceId>> face_cells_map;
- for (unsigned int j=0; j<m_number_of_cells; ++j) {
- const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
- for (unsigned short r=0; r<cell_nb_nodes; ++r) {
- unsigned int node0_id = m_cell_nodes(j,r);
- unsigned int node1_id = m_cell_nodes(j,(r+1)%cell_nb_nodes);
- if (node1_id<node0_id) {
- std::swap(node0_id, node1_id);
- }
- face_cells_map[Face(node0_id, node1_id)].push_back(std::make_pair(j, r));
- }
- }
-
- m_number_of_faces = face_cells_map.size();
- Kokkos::View<unsigned short*> face_nb_nodes("face_nb_nodes", m_number_of_faces);
- Kokkos::parallel_for(m_number_of_faces, KOKKOS_LAMBDA(const unsigned int& l) {
- face_nb_nodes[l] = 2;
- });
- m_face_nb_nodes = face_nb_nodes;
-
- Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, 2);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const Face& face = face_cells_vector.first;
- face_nodes(l,0) = face.m_node0_id;
- face_nodes(l,1) = face.m_node1_id;
- ++l;
- }
- }
- m_face_nodes = face_nodes;
-
- Kokkos::View<unsigned short*> face_nb_cells("face_nb_cells", m_number_of_faces);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- face_nb_cells[l] = cells_vector.size();
- ++l;
- }
- }
- m_face_nb_cells = face_nb_cells;
-
- Kokkos::View<unsigned int**> face_cells("face_cells", face_cells_map.size(), 2);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned int cell_number = cells_vector[lj].first;
- face_cells(l,lj) = cell_number;
- }
- ++l;
- }
- }
- m_face_cells = face_cells;
-
- // In 2d cell_nb_faces = cell_nb_node
- m_cell_nb_faces = m_cell_nb_nodes;
- Kokkos::View<unsigned int**> cell_faces("cell_faces", m_number_of_faces, m_max_nb_node_per_cell);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned int cell_number = cells_vector[lj].first;
- unsigned short cell_local_face = cells_vector[lj].second;
- cell_faces(cell_number,cell_local_face) = l;
- }
- ++l;
- }
- }
- m_cell_faces = cell_faces;
-
- Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
- m_number_of_faces, m_max_nb_node_per_cell);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned short cell_local_face = cells_vector[lj].second;
- face_cell_local_face(l,lj) = cell_local_face;
- }
- ++l;
- }
- }
- m_face_cell_local_face = face_cell_local_face;
- }
-
- public:
- void addRefFaceList(const RefFaceList& ref_face_list)
- {
- m_ref_face_list.push_back(ref_face_list);
- }
-
- size_t numberOfRefFaceList() const
- {
- return m_ref_face_list.size();
- }
-
- const RefFaceList& refFaceList(const size_t& i) const
- {
- return m_ref_face_list[i];
- }
-
- void addRefNodeList(const RefNodeList& ref_node_list)
- {
- m_ref_node_list.push_back(ref_node_list);
- }
-
- size_t numberOfRefNodeList() const
- {
- return m_ref_node_list.size();
- }
-
- const RefNodeList& refNodeList(const size_t& i) const
- {
- return m_ref_node_list[i];
- }
-
- const size_t& numberOfNodes() const
- {
- return m_number_of_nodes;
- }
-
- const size_t& numberOfFaces() const
- {
- return m_number_of_faces;
- }
-
- const size_t& numberOfCells() const
- {
- return m_number_of_cells;
- }
-
- const size_t& maxNbNodePerCell() const
- {
- return m_max_nb_node_per_cell;
- }
-
- const Kokkos::View<const unsigned int**> cellNodes() const
- {
- return m_cell_nodes;
- }
-
- const Kokkos::View<const unsigned int**> cellFaces() const
- {
- return m_cell_faces;
- }
-
- const Kokkos::View<const unsigned short*> nodeNbCells() const
- {
- return m_node_nb_cells;
- }
-
- const Kokkos::View<const unsigned short*> cellNbNodes() const
- {
- return m_cell_nb_nodes;
- }
-
- const Kokkos::View<const double*> invCellNbNodes() const
- {
- return m_inv_cell_nb_nodes;
- }
-
- const Kokkos::View<const unsigned short*> cellNbFaces() const
- {
- return m_cell_nb_faces;
- }
-
- const Kokkos::View<const unsigned short*> faceNbCells() const
- {
- return m_face_nb_cells;
- }
-
- const Kokkos::View<const unsigned short*> faceNbNodes() const
- {
- return m_face_nb_nodes;
- }
-
- const Kokkos::View<const unsigned int**> nodeCells() const
- {
- return m_node_cells;
- }
-
- const Kokkos::View<const unsigned int**> faceCells() const
- {
- return m_face_cells;
- }
-
- const Kokkos::View<const unsigned int**> faceNodes() const
- {
- return m_face_nodes;
- }
-
- const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
- {
- return m_node_cell_local_node;
- }
-
- const Kokkos::View<const unsigned short**> faceCellLocalFace() const
- {
- return m_face_cell_local_face;
- }
-
- unsigned int getFaceNumber(const unsigned int node0_id,
- const unsigned int node1_id) const
- {
-#warning rewrite
- const unsigned int n0_id = std::min(node0_id, node1_id);
- const unsigned int n1_id = std::max(node0_id, node1_id);
- // worst code ever
- for (unsigned int l=0; l<m_number_of_faces; ++l) {
- if ((m_face_nodes(l,0) == n0_id) and (m_face_nodes(l,1) == n1_id)) {
- return l;
- }
- }
- std::cerr << "Face not found!\n";
- std::exit(0);
- return 0;
- }
-
- Connectivity2D(const Connectivity2D&) = delete;
-
- Connectivity2D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
- const Kokkos::View<const unsigned int**> cell_nodes)
- : m_number_of_cells (cell_nodes.extent(0)),
- m_cell_nb_nodes(cell_nb_nodes),
- m_cell_nodes (cell_nodes)
- {
- {
- Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){
- inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
- });
- m_inv_cell_nb_nodes = inv_cell_nb_nodes;
- }
- Assert(m_number_of_cells>0);
-
- ConnectivityUtils utils;
- utils.computeNodeCellConnectivity(m_cell_nodes,
- m_cell_nb_nodes,
- m_number_of_cells,
- m_max_nb_node_per_cell,
- m_number_of_nodes,
- m_node_nb_cells,
- m_node_cells,
- m_node_cell_local_node);
-
- this->_computeFaceCellConnectivities();
- }
-
- ~Connectivity2D()
- {
- ;
- }
-};
-
-#endif // CONNECTIVITY_2D_HPP
diff --git a/src/mesh/Connectivity3D.hpp b/src/mesh/Connectivity3D.hpp
deleted file mode 100644
index ce982d2765dd82736c3211238017529ce997d92f..0000000000000000000000000000000000000000
--- a/src/mesh/Connectivity3D.hpp
+++ /dev/null
@@ -1,566 +0,0 @@
-#ifndef CONNECTIVITY_3D_HPP
-#define CONNECTIVITY_3D_HPP
-
-#include <Kokkos_Core.hpp>
-#include <PastisAssert.hpp>
-#include <TinyVector.hpp>
-
-#include <ConnectivityUtils.hpp>
-#include <vector>
-#include <map>
-#include <unordered_map>
-#include <algorithm>
-
-#include <RefId.hpp>
-#include <RefNodeList.hpp>
-#include <RefFaceList.hpp>
-
-#include <tuple>
-
-class Connectivity3D
-{
-public:
- static constexpr size_t dimension = 3;
-
-private:
- std::vector<RefFaceList> m_ref_face_list;
- std::vector<RefNodeList> m_ref_node_list;
-
- const size_t m_number_of_cells;
- size_t m_number_of_faces;
- size_t m_number_of_nodes;
-
- const Kokkos::View<const unsigned short*> m_cell_nb_nodes;
- const Kokkos::View<const unsigned int**> m_cell_nodes;
- Kokkos::View<double*> m_inv_cell_nb_nodes;
-
- Kokkos::View<const unsigned short*> m_cell_nb_faces;
- Kokkos::View<const unsigned int**> m_cell_faces;
- Kokkos::View<const bool**> m_cell_faces_is_reversed;
-
- Kokkos::View<const unsigned short*> m_node_nb_cells;
- Kokkos::View<const unsigned int**> m_node_cells;
- Kokkos::View<const unsigned short**> m_node_cell_local_node;
-
- Kokkos::View<const unsigned short*> m_face_nb_cells;
- Kokkos::View<const unsigned int**> m_face_cells;
- Kokkos::View<const unsigned short**> m_face_cell_local_face;
-
- Kokkos::View<const unsigned short*> m_face_nb_nodes;
- Kokkos::View<const unsigned int**> m_face_nodes;
- // Kokkos::View<unsigned short**> m_face_node_local_face;
-
- Kokkos::View<const unsigned short*> m_node_nb_faces;
- Kokkos::View<const unsigned int**> m_node_faces;
- // Kokkos::View<unsigned short**> m_node_face_local_node;
-
- size_t m_max_nb_node_per_cell;
- size_t m_max_nb_face_per_cell;
- size_t m_max_nb_node_per_face;
-
- class Face
- {
- public:
- friend struct Hash;
- struct Hash {
- size_t operator()(const Face& f) const
- {
- size_t hash = 0;
- for (size_t i=0; i<f.m_node_id_list.size(); ++i) {
- hash ^= std::hash<unsigned int>()(f.m_node_id_list[i]) >> i;
- }
- return hash;
- }
- };
- private:
- bool m_reversed;
- std::vector<unsigned int> m_node_id_list;
-
- public:
- friend std::ostream& operator<<(std::ostream& os, const Face& f)
- {
- for (auto id : f.m_node_id_list) {
- std::cout << id << ' ';
- }
- return os;
- }
-
- KOKKOS_INLINE_FUNCTION
- const bool& reversed() const
- {
- return m_reversed;
- }
-
- KOKKOS_INLINE_FUNCTION
- const std::vector<unsigned int>& nodeIdList() const
- {
- return m_node_id_list;
- }
-
- KOKKOS_INLINE_FUNCTION
- std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
- {
- const auto min_id = std::min_element(node_list.begin(), node_list.end());
- const int shift = std::distance(node_list.begin(), min_id);
-
- std::vector<unsigned int> rotated_node_list(node_list.size());
- if (node_list[(shift+1)%node_list.size()] > node_list[(shift+node_list.size()-1)%node_list.size()]) {
- for (size_t i=0; i<node_list.size(); ++i) {
- rotated_node_list[i] = node_list[(shift+node_list.size()-i)%node_list.size()];
- m_reversed = true;
- }
- } else {
- for (size_t i=0; i<node_list.size(); ++i) {
- rotated_node_list[i] = node_list[(shift+i)%node_list.size()];
- }
- }
-
- return rotated_node_list;
- }
-
- bool operator==(const Face& f) const
- {
- if (m_node_id_list.size() == f.nodeIdList().size()) {
- for (size_t j=0; j<m_node_id_list.size(); ++j) {
- if (m_node_id_list[j] != f.nodeIdList()[j]) {
- return false;
- }
- }
- return true;
- }
- return false;
- }
-
- KOKKOS_INLINE_FUNCTION
- bool operator<(const Face& f) const
- {
- const size_t min_nb_nodes = std::min(f.m_node_id_list.size(), m_node_id_list.size());
- for (size_t i=0; i<min_nb_nodes; ++i) {
- if (m_node_id_list[i] < f.m_node_id_list[i]) return true;
- if (m_node_id_list[i] != f.m_node_id_list[i]) return false;
- }
- return m_node_id_list.size() < f.m_node_id_list.size();
- }
-
- KOKKOS_INLINE_FUNCTION
- Face& operator=(const Face&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face& operator=(Face&&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(const Face&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(Face&&) = default;
-
- KOKKOS_INLINE_FUNCTION
- Face(const std::vector<unsigned int>& given_node_id_list)
- : m_reversed(false),
- m_node_id_list(_sort(given_node_id_list))
- {
- ;
- }
-
- KOKKOS_INLINE_FUNCTION
- Face() = delete;
-
- KOKKOS_INLINE_FUNCTION
- ~Face() = default;
- };
-
- std::unordered_map<Face,unsigned int, Face::Hash> m_face_number_map;
-
- void _computeFaceCellConnectivities()
- {
- Kokkos::View<unsigned short*> cell_nb_faces("cell_nb_faces", m_number_of_cells);
-
- typedef std::tuple<unsigned int, unsigned short, bool> CellFaceInfo;
- std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
- for (unsigned int j=0; j<m_number_of_cells; ++j) {
- const unsigned short cell_nb_nodes = m_cell_nb_nodes[j];
- switch (cell_nb_nodes) {
- case 4: { // tetrahedron
- cell_nb_faces[j] = 4;
- // face 0
- Face f0({m_cell_nodes(j,1),
- m_cell_nodes(j,2),
- m_cell_nodes(j,3)});
- face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({m_cell_nodes(j,0),
- m_cell_nodes(j,3),
- m_cell_nodes(j,2)});
- face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({m_cell_nodes(j,0),
- m_cell_nodes(j,1),
- m_cell_nodes(j,3)});
- face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed()));
-
- // face 3
- Face f3({m_cell_nodes(j,0),
- m_cell_nodes(j,2),
- m_cell_nodes(j,1)});
- face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed()));
- break;
- }
- case 8: { // hexahedron
- // face 0
- Face f0({m_cell_nodes(j,3),
- m_cell_nodes(j,2),
- m_cell_nodes(j,1),
- m_cell_nodes(j,0)});
- face_cells_map[f0].push_back(std::make_tuple(j, 0, f0.reversed()));
-
- // face 1
- Face f1({m_cell_nodes(j,4),
- m_cell_nodes(j,5),
- m_cell_nodes(j,6),
- m_cell_nodes(j,7)});
- face_cells_map[f1].push_back(std::make_tuple(j, 1, f1.reversed()));
-
- // face 2
- Face f2({m_cell_nodes(j,0),
- m_cell_nodes(j,4),
- m_cell_nodes(j,7),
- m_cell_nodes(j,3)});
- face_cells_map[f2].push_back(std::make_tuple(j, 2, f2.reversed()));
-
- // face 3
- Face f3({m_cell_nodes(j,1),
- m_cell_nodes(j,2),
- m_cell_nodes(j,6),
- m_cell_nodes(j,5)});
- face_cells_map[f3].push_back(std::make_tuple(j, 3, f3.reversed()));
-
- // face 4
- Face f4({m_cell_nodes(j,0),
- m_cell_nodes(j,1),
- m_cell_nodes(j,5),
- m_cell_nodes(j,4)});
- face_cells_map[f4].push_back(std::make_tuple(j, 4, f4.reversed()));
-
- // face 5
- Face f5({m_cell_nodes(j,3),
- m_cell_nodes(j,7),
- m_cell_nodes(j,6),
- m_cell_nodes(j,2)});
- face_cells_map[f5].push_back(std::make_tuple(j, 5, f5.reversed()));
-
- cell_nb_faces[j] = 6;
- break;
- }
- default: {
- std::cerr << "unexpected cell type!\n";
- std::exit(0);
- }
- }
- }
- std::cerr << __FILE__ << ':' << __LINE__ << ':'
- << rang::fg::red << " only works for hexa" << rang::fg::reset << '\n';
-
- m_cell_nb_faces = cell_nb_faces;
-
- m_number_of_faces = face_cells_map.size();
-
- std::cerr << __FILE__ << ':' << __LINE__ << ':'
- << rang::fg::red << " m_max_nb_node_per_face forced to 4" << rang::fg::reset << '\n';
- Kokkos::View<unsigned short*> face_nb_nodes("face_nb_nodes", m_number_of_faces);
- m_max_nb_node_per_face = 4;
- Kokkos::View<unsigned int**> face_nodes("face_nodes", m_number_of_faces, m_max_nb_node_per_face);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const Face& face = face_cells_vector.first;
- for(size_t r=0; r<face.nodeIdList().size(); ++r) {
- face_nodes(l,r) = face.nodeIdList()[r];
- }
- face_nb_nodes[l] = face.nodeIdList().size();
- m_face_number_map[face] = l;
- ++l;
- }
- }
- m_face_nodes = face_nodes;
- m_face_nb_nodes = face_nb_nodes;
-
- Kokkos::View<unsigned short*> face_nb_cells("face_nb_cells", m_number_of_faces);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- face_nb_cells[l] = cells_vector.size();
- ++l;
- }
- }
- m_face_nb_cells = face_nb_cells;
-
-#warning check that the number of cell per faces is <=2
-
- Kokkos::View<unsigned int**> face_cells("face_cells", face_cells_map.size(), 2);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- const unsigned int cell_number = std::get<0>(cells_vector[lj]);
- face_cells(l,lj) = cell_number;
- }
- ++l;
- }
- }
- m_face_cells = face_cells;
-
- std::cerr << __FILE__ << ':' << __LINE__ << ':'
- << rang::fg::red << " m_max_nb_face_per_cell forced to 6" << rang::fg::reset << '\n';
- m_max_nb_face_per_cell = 6;
- Kokkos::View<unsigned int**> cell_faces("cell_faces", m_number_of_cells, m_max_nb_face_per_cell);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
- cell_faces(cell_number,cell_local_face) = l;
- }
- ++l;
- }
- }
- m_cell_faces = cell_faces;
-
- Kokkos::View<bool**> cell_faces_is_reversed("cell_faces_is_reversed", m_number_of_cells, m_max_nb_face_per_cell);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
- cell_faces_is_reversed(cell_number,cell_local_face) = reversed;
- }
- ++l;
- }
- }
- m_cell_faces_is_reversed = cell_faces_is_reversed;
-
- Kokkos::View<unsigned short**> face_cell_local_face("face_cell_local_face",
- m_number_of_faces, m_max_nb_node_per_cell);
- {
- int l=0;
- for (const auto& face_cells_vector : face_cells_map) {
- const auto& cells_vector = face_cells_vector.second;
- for (unsigned short lj=0; lj<cells_vector.size(); ++lj) {
- unsigned short cell_local_face = std::get<1>(cells_vector[lj]);
- face_cell_local_face(l,lj) = cell_local_face;
- }
- ++l;
- }
- }
- m_face_cell_local_face = face_cell_local_face;
-
- std::unordered_map<unsigned int, std::vector<unsigned int>> node_faces_map;
- for (size_t l=0; l<m_number_of_faces; ++l) {
- for (size_t lr=0; lr<face_nb_nodes(l); ++lr) {
- const unsigned int r = face_nodes(l, lr);
- node_faces_map[r].push_back(l);
- }
- }
- Kokkos::View<unsigned short*> node_nb_faces("node_nb_faces", m_number_of_nodes);
- size_t max_nb_face_per_node = 0;
- for (auto node_faces : node_faces_map) {
- max_nb_face_per_node = std::max(node_faces.second.size(), max_nb_face_per_node);
- node_nb_faces[node_faces.first] = node_faces.second.size();
- }
- m_node_nb_faces = node_nb_faces;
-
- Kokkos::View<unsigned int**> node_faces("node_faces", m_number_of_nodes, max_nb_face_per_node);
- for (auto node_faces_vector : node_faces_map) {
- const unsigned int r = node_faces_vector.first;
- const std::vector<unsigned int>& faces_vector = node_faces_vector.second;
- for (size_t l=0; l < faces_vector.size(); ++l) {
- node_faces(r, l) = faces_vector[l];
- }
- }
- m_node_faces = node_faces;
- }
-
- public:
- void addRefFaceList(const RefFaceList& ref_face_list)
- {
- m_ref_face_list.push_back(ref_face_list);
- }
-
- size_t numberOfRefFaceList() const
- {
- return m_ref_face_list.size();
- }
-
- const RefFaceList& refFaceList(const size_t& i) const
- {
- return m_ref_face_list[i];
- }
-
- void addRefNodeList(const RefNodeList& ref_node_list)
- {
- m_ref_node_list.push_back(ref_node_list);
- }
-
- size_t numberOfRefNodeList() const
- {
- return m_ref_node_list.size();
- }
-
- const RefNodeList& refNodeList(const size_t& i) const
- {
- return m_ref_node_list[i];
- }
-
- const size_t& numberOfNodes() const
- {
- return m_number_of_nodes;
- }
-
- const size_t& numberOfFaces() const
- {
- return m_number_of_faces;
- }
-
- const size_t& numberOfCells() const
- {
- return m_number_of_cells;
- }
-
- const size_t& maxNbFacePerCell() const
- {
- return m_max_nb_face_per_cell;
- }
-
- const size_t& maxNbNodePerCell() const
- {
- return m_max_nb_node_per_cell;
- }
-
- const size_t& maxNbNodePerFace() const
- {
- return m_max_nb_node_per_face;
- }
-
- const Kokkos::View<const unsigned int**> cellNodes() const
- {
- return m_cell_nodes;
- }
-
- const Kokkos::View<const unsigned int**> cellFaces() const
- {
- return m_cell_faces;
- }
-
- const Kokkos::View<const bool**> cellFacesIsReversed() const
- {
- return m_cell_faces_is_reversed;
- }
-
- const Kokkos::View<const unsigned short*> nodeNbCells() const
- {
- return m_node_nb_cells;
- }
-
- const Kokkos::View<const unsigned short*> cellNbNodes() const
- {
- return m_cell_nb_nodes;
- }
-
- const Kokkos::View<const double*> invCellNbNodes() const
- {
- return m_inv_cell_nb_nodes;
- }
-
- const Kokkos::View<const unsigned short*> cellNbFaces() const
- {
- return m_cell_nb_faces;
- }
-
- const Kokkos::View<const unsigned short*> faceNbCells() const
- {
- return m_face_nb_cells;
- }
-
- const Kokkos::View<const unsigned short*> faceNbNodes() const
- {
- return m_face_nb_nodes;
- }
-
- const Kokkos::View<const unsigned int**> nodeCells() const
- {
- return m_node_cells;
- }
-
- const Kokkos::View<const unsigned int**> faceCells() const
- {
- return m_face_cells;
- }
-
- const Kokkos::View<const unsigned int**> faceNodes() const
- {
- return m_face_nodes;
- }
-
- const Kokkos::View<const unsigned short**> nodeCellLocalNode() const
- {
- return m_node_cell_local_node;
- }
-
- const Kokkos::View<const unsigned short**> faceCellLocalFace() const
- {
- return m_face_cell_local_face;
- }
-
- unsigned int getFaceNumber(const std::vector<unsigned int>& face_nodes) const
- {
- const Face face(face_nodes);
- auto i_face = m_face_number_map.find(face);
- if (i_face == m_face_number_map.end()) {
- std::cerr << "Face " << face << " not found!\n";
- std::exit(0);
- }
- return i_face->second;
- }
-
- Connectivity3D(const Connectivity3D&) = delete;
-
- Connectivity3D(const Kokkos::View<const unsigned short*> cell_nb_nodes,
- const Kokkos::View<const unsigned int**> cell_nodes)
- : m_number_of_cells (cell_nodes.extent(0)),
- m_cell_nb_nodes(cell_nb_nodes),
- m_cell_nodes (cell_nodes)
- {
- {
- Kokkos::View<double*> inv_cell_nb_nodes("inv_cell_nb_nodes", m_number_of_cells);
- Kokkos::parallel_for(m_number_of_cells, KOKKOS_LAMBDA(const int&j){
- inv_cell_nb_nodes[j] = 1./m_cell_nb_nodes[j];
- });
- m_inv_cell_nb_nodes = inv_cell_nb_nodes;
- }
- Assert(m_number_of_cells>0);
-
- ConnectivityUtils utils;
- utils.computeNodeCellConnectivity(m_cell_nodes,
- m_cell_nb_nodes,
- m_number_of_cells,
- m_max_nb_node_per_cell,
- m_number_of_nodes,
- m_node_nb_cells,
- m_node_cells,
- m_node_cell_local_node);
-
- this->_computeFaceCellConnectivities();
- }
-
- ~Connectivity3D()
- {
- ;
- }
-};
-
-#endif // CONNECTIVITY_3D_HPP
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0c87ccb6d64f840cd450eb5ded67c6a337a5cf25
--- /dev/null
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -0,0 +1,158 @@
+#include <Connectivity.hpp>
+
+#include <ConnectivityComputer.hpp>
+#include <map>
+
+
+template <typename ConnectivityType>
+KOKKOS_INLINE_FUNCTION
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const ConnectivityType& connectivity,
+ const ItemType& item_type,
+ const ItemType& child_item_type) const
+{
+ 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);
+
+ std::cout << "computing connectivity "
+ << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
+
+ item_to_child_item_matrix
+ = this->_computeInverse(child_to_item_matrix);
+ } else {
+ std::cerr << "unable to compute connectivity "
+ << itemName(item_type) << " -> " << itemName(child_item_type) << '\n';
+ std::exit(0);
+ }
+
+ return item_to_child_item_matrix;
+}
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity1D&, const ItemType&, const ItemType&) const;
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity2D&, const ItemType&, const ItemType&) const;
+
+
+template
+ConnectivityMatrix ConnectivityComputer::
+computeConnectivityMatrix(const Connectivity3D&, const ItemType&, const ItemType&) const;
+
+ConnectivityMatrix
+ConnectivityComputer::
+_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
+{
+ std::map<unsigned int, std::vector<unsigned int>> child_to_item_vector_map;
+ const size_t& number_of_items = item_to_child_matrix.numRows();
+
+ for (unsigned int j=0; j<number_of_items; ++j) {
+ const auto& item_to_child = item_to_child_matrix.rowConst(j);
+ for (unsigned int r=0; r<item_to_child.length; ++r) {
+ child_to_item_vector_map[item_to_child(r)].push_back(j);
+ }
+ }
+
+ {
+ size_t i=0;
+ for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
+ if (child_item_id != i) {
+ std::cerr << "sparse item numerotation NIY\n";
+ std::exit(0);
+ }
+ ++i;
+ }
+ }
+
+ std::vector<std::vector<unsigned int>> child_to_items_vector(child_to_item_vector_map.size());
+ for (const auto& [child_item_id, item_vector] : child_to_item_vector_map) {
+ child_to_items_vector[child_item_id] = item_vector;
+ }
+
+ return ConnectivityMatrix(child_to_items_vector);
+}
+
+template <ItemType item_type,
+ ItemType child_item_type,
+ typename ConnectivityType>
+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);
+ const ConnectivityMatrix& item_to_child_items_matrix
+ = connectivity.getMatrix(item_type, child_item_type);
+
+ SubItemValuePerItem<unsigned short, child_item_type, item_type> item_number_in_child_item(connectivity);
+ for (unsigned int r=0; r<item_to_child_items_matrix.numRows(); ++r) {
+ const auto& item_to_child_items = item_to_child_items_matrix.rowConst(r);
+ for (unsigned short J=0; J<item_to_child_items.length; ++J) {
+ const unsigned int j = item_to_child_items(J);
+ const auto& child_item_to_items = child_item_to_items_matrix.rowConst(j);
+
+ for (unsigned int R=0; R<child_item_to_items.length; ++R) {
+ if (child_item_to_items(R) == r) {
+ item_number_in_child_item(r,J) = R;
+ break;
+ }
+ }
+ }
+ }
+
+ return item_number_in_child_item;
+}
+
+// 1D
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::node,
+ ItemType::cell>(const Connectivity1D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::face,
+ ItemType::cell>(const Connectivity1D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::cell,
+ ItemType::node>(const Connectivity1D&) const;
+
+// 2D
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::node,
+ ItemType::cell>(const Connectivity2D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::face,
+ ItemType::cell>(const Connectivity2D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::cell,
+ ItemType::node>(const Connectivity2D&) const;
+
+// 3D
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::node>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::node,
+ ItemType::cell>(const Connectivity3D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::cell, ItemType::face>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::face,
+ ItemType::cell>(const Connectivity3D&) const;
+
+template SubItemValuePerItem<const unsigned short, ItemType::node, ItemType::cell>
+ConnectivityComputer::
+computeLocalItemNumberInChildItem<ItemType::cell,
+ ItemType::node>(const Connectivity3D&) const;
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..8b4fdfc3ab20bf2d6290d234a13b833dc6451cb6
--- /dev/null
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -0,0 +1,32 @@
+#ifndef CONNECTIVITY_COMPUTER_HPP
+#define CONNECTIVITY_COMPUTER_HPP
+
+#include <ConnectivityMatrix.hpp>
+#include <SubItemValuePerItem.hpp>
+
+class ConnectivityComputer
+{
+ private:
+ ConnectivityMatrix
+ _computeInverse(const ConnectivityMatrix& item_to_child_matrix) const;
+
+ public:
+ template <typename ConnectivityType>
+ ConnectivityMatrix
+ computeConnectivityMatrix(const ConnectivityType& connectivity,
+ const ItemType& item_type,
+ const ItemType& child_item_type) const;
+
+
+ template <ItemType item_type,
+ ItemType child_item_type,
+ typename ConnectivityType>
+ SubItemValuePerItem<const unsigned short, child_item_type, item_type>
+ computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const;
+
+ ConnectivityComputer(const ConnectivityComputer&) = default;
+ ConnectivityComputer() = default;
+ ~ConnectivityComputer() = default;
+};
+
+#endif // CONNECTIVITY_COMPUTER_HPP
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..356ee6da56d474836180cc867318c2d4f912549a
--- /dev/null
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -0,0 +1,74 @@
+#ifndef CONNECTIVITY_MATRIX_HPP
+#define CONNECTIVITY_MATRIX_HPP
+
+#include <Kokkos_Core.hpp>
+#include <Kokkos_StaticCrsGraph.hpp>
+
+class ConnectivityMatrix
+{
+ private:
+ typedef Kokkos::StaticCrsGraph<unsigned int, Kokkos::HostSpace> HostMatrix;
+ HostMatrix m_host_matrix;
+
+ bool m_is_built{false};
+
+ public:
+ typedef HostMatrix::row_map_type HostRowType;
+
+ const bool& isBuilt() const
+ {
+ return m_is_built;
+ }
+
+ const auto& entries() const
+ {
+ return m_host_matrix.entries;
+ }
+
+ const auto& rowsMap() const
+ {
+ return m_host_matrix.row_map;
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto numEntries() const
+ {
+ return m_host_matrix.entries.extent(0);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto numRows() const
+ {
+ return m_host_matrix.numRows();
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto rowConst(const size_t& j) const
+ {
+ return m_host_matrix.rowConst(j);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const auto rowMap(const size_t& j) const
+ {
+ return m_host_matrix.row_map[j];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer)
+ : m_host_matrix{Kokkos::create_staticcrsgraph<HostMatrix>("connectivity_matrix", initializer)},
+ m_is_built{true}
+ {
+ ;
+ }
+
+ ConnectivityMatrix& operator=(const ConnectivityMatrix&) = default;
+ ConnectivityMatrix& operator=(ConnectivityMatrix&&) = default;
+
+ ConnectivityMatrix() = default;
+ ConnectivityMatrix(const ConnectivityMatrix&) = default;
+ ConnectivityMatrix(ConnectivityMatrix&&) = default;
+ ~ConnectivityMatrix() = default;
+};
+
+#endif // CONNECTIVITY_MATRIX_HPP
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
deleted file mode 100644
index 3fe9e6f6a11032e752ce5bdcd9d76e7a9c46e1b0..0000000000000000000000000000000000000000
--- a/src/mesh/ConnectivityUtils.hpp
+++ /dev/null
@@ -1,84 +0,0 @@
-#ifndef CONNECTIVITY_UTILS_HPP
-#define CONNECTIVITY_UTILS_HPP
-
-#include <map>
-#include <Kokkos_Core.hpp>
-
-class ConnectivityUtils
-{
- public:
- void computeNodeCellConnectivity(const Kokkos::View<const unsigned int**> cell_nodes,
- const Kokkos::View<const unsigned short*> cell_nb_nodes,
- const size_t& number_of_cells,
- size_t& max_nb_node_per_cell,
- size_t& number_of_nodes,
- Kokkos::View<const unsigned short*>& node_nb_cells,
- Kokkos::View<const unsigned int**>& node_cells,
- Kokkos::View<const unsigned short**>& node_cell_local_node)
- {
- std::map<unsigned int, std::vector<unsigned int>> node_cells_map;
- using namespace Kokkos::Experimental;
- Kokkos::parallel_reduce(number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
- const size_t n = cell_nb_nodes[j];
- if (n > nb_max) nb_max = n;
- }, Kokkos::Max<size_t>(max_nb_node_per_cell));
-
- for (unsigned int j=0; j<number_of_cells; ++j) {
- for (unsigned int r=0; r<cell_nb_nodes[j]; ++r) {
- node_cells_map[cell_nodes(j,r)].push_back(j);
- }
- }
-
- {
- size_t i=0;
- for (const auto& i_cell_vector : node_cells_map) {
- const auto& [node_id, cell_vector] = i_cell_vector;
- if (node_id != i) {
- std::cerr << "sparse node numerotation NIY\n";
- std::exit(0);
- }
- ++i;
- }
- }
-
- number_of_nodes = node_cells_map.size();
-
- Kokkos::View<unsigned short*> built_node_nb_cells("node_nb_cells", node_cells_map.size());
- size_t max_node_cells = 0;
- for (const auto& i_cell_vector : node_cells_map) {
- const auto& [i, cells_vector] = i_cell_vector;
- const size_t nb_cells = cells_vector.size();
- built_node_nb_cells[i] = nb_cells;
- if (nb_cells > max_node_cells) {
- max_node_cells = nb_cells;
- }
- }
- node_nb_cells = built_node_nb_cells;
-
- Kokkos::View<unsigned int**> built_node_cells("node_cells", node_cells_map.size(), max_node_cells);
- for (const auto& i_cell_vector : node_cells_map) {
- const auto& [i, cells_vector] = i_cell_vector;
- for (size_t j=0; j<cells_vector.size(); ++j) {
- built_node_cells(i,j) = cells_vector[j];
- }
- }
- node_cells = built_node_cells;
-
- Kokkos::View<unsigned short**> built_node_cell_local_node("node_cell_local_node",
- node_cells_map.size(), max_node_cells);
- Kokkos::parallel_for(number_of_nodes, KOKKOS_LAMBDA(const unsigned int& r){
- for (unsigned short J=0; J<node_nb_cells[r]; ++J) {
- const unsigned int j = node_cells(r,J);
- for (unsigned int R=0; R<cell_nb_nodes[j]; ++R) {
- if (cell_nodes(j,R) == r) {
- built_node_cell_local_node(r,J)=R;
- break;
- }
- }
- }
- });
- node_cell_local_node = built_node_cell_local_node;
- }
-};
-
-#endif // CONNECTIVITY_UTILS_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index c8c2e09e7a28bf484cac0840b41d964f87fc994a..f16ab2fa02e06aad3715eeb7c84028476c7565ae 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -5,9 +5,8 @@
#include <set>
#include <rang.hpp>
-#include <Connectivity1D.hpp>
-#include <Connectivity2D.hpp>
-#include <Connectivity3D.hpp>
+#include <CellType.hpp>
+#include <Connectivity.hpp>
#include <Mesh.hpp>
@@ -52,7 +51,7 @@ class ErrorHandler
* Prints the error message
*
*/
- virtual void writeErrorMessage();
+ virtual void writeErrorMessage() const;
/**
* The copy constructor
@@ -90,7 +89,7 @@ class ErrorHandler
;
}
};
-void ErrorHandler::writeErrorMessage()
+void ErrorHandler::writeErrorMessage() const
{
switch(__type) {
case asked: {
@@ -293,7 +292,7 @@ GmshReader::GmshReader(const std::string& filename)
this->__proceedData();
// this->__createMesh();
}
- catch(ErrorHandler e) {
+ catch(const ErrorHandler& e) {
e.writeErrorMessage();
std::exit(0);
}
@@ -789,39 +788,36 @@ GmshReader::__proceedData()
dimension3_mask[12]=1;
dimension3_mask[13]=1;
- std::cout << "- dimension 0 entities : " << (dimension0_mask, elementNumber) << '\n';
- std::cout << "- dimension 1 entities : " << (dimension1_mask, elementNumber) << '\n';
- std::cout << "- dimension 2 entities : " << (dimension2_mask, elementNumber) << '\n';
- std::cout << "- dimension 3 entities : " << (dimension3_mask, elementNumber) << '\n';
+ std::cout << "- dimension 0 entities: " << (dimension0_mask, elementNumber) << '\n';
+ std::cout << "- dimension 1 entities: " << (dimension1_mask, elementNumber) << '\n';
+ std::cout << "- dimension 2 entities: " << (dimension2_mask, elementNumber) << '\n';
+ std::cout << "- dimension 3 entities: " << (dimension3_mask, elementNumber) << '\n';
if ((dimension3_mask, elementNumber)>0) {
const size_t nb_cells = (dimension3_mask, elementNumber);
- size_t max_nb_node_per_cell=4;
- if (elementNumber[4] > 0) {
- max_nb_node_per_cell = 8;
- }
- const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+
+ std::vector<CellType> cell_type_vector(nb_cells);
+
+ std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
const size_t nb_tetrahedra = __tetrahedra.extent(0);
for (size_t j=0; j<nb_tetrahedra; ++j) {
+ cell_by_node_vector[j].resize(4);
for (int r=0; r<4; ++r) {
- cell_nodes(j,r) = __tetrahedra[j][r];
+ cell_by_node_vector[j][r] = __tetrahedra[j][r];
}
+ cell_type_vector[j] = CellType::Tetrahedron;
}
const size_t nb_hexahedra = __hexahedra.extent(0);
for (size_t j=0; j<nb_hexahedra; ++j) {
- const size_t jh = j+nb_tetrahedra;
+ const size_t jh = nb_tetrahedra+j;
+ cell_by_node_vector[jh].resize(8);
for (int r=0; r<8; ++r) {
- cell_nodes(jh,r) = __hexahedra[j][r];
+ cell_by_node_vector[jh][r] = __hexahedra[j][r];
}
- }
- const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
- for (size_t j=0; j<nb_tetrahedra; ++j) {
- cell_nb_nodes[j] = 4;
- }
- for (size_t j=nb_tetrahedra; j<nb_tetrahedra+nb_hexahedra; ++j) {
- cell_nb_nodes[j] = 8;
+ cell_type_vector[jh] = CellType::Hexahedron;
}
- std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_nb_nodes, cell_nodes));
+ std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector,
+ cell_type_vector));
Connectivity3D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -847,16 +843,6 @@ GmshReader::__proceedData()
connectivity.addRefFaceList(RefFaceList(physical_ref_id.refId(), face_list));
}
- // for (size_t j=0; j<nb_cells; ++j) {
- // std::cout << std::setw(5) << j << ": ";
- // for (size_t r=0; r<cell_nb_nodes[j]; ++r) {
- // std::cout << ' ' << cell_nodes(j,r);
- // }
- // std::cout << '\n';
- // }
-
- // std::cerr << "*** using a 3d mesh (NIY)\n";
- // std::exit(0);
typedef Mesh<Connectivity3D> MeshType;
typedef TinyVector<3, double> Rd;
@@ -870,39 +856,36 @@ GmshReader::__proceedData()
} else if ((dimension2_mask, elementNumber)>0) {
const size_t nb_cells = (dimension2_mask, elementNumber);
- size_t max_nb_node_per_cell=3;
- if (elementNumber[2] > 0) {
- max_nb_node_per_cell = 4;
- }
- const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
+
+ std::vector<CellType> cell_type_vector(nb_cells);
+
+ std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
const size_t nb_triangles = __triangles.extent(0);
for (size_t j=0; j<nb_triangles; ++j) {
- cell_nodes(j,0) = __triangles[j][0];
- cell_nodes(j,1) = __triangles[j][1];
- cell_nodes(j,2) = __triangles[j][2];
+ cell_by_node_vector[j].resize(3);
+ for (int r=0; r<3; ++r) {
+ cell_by_node_vector[j][r] = __triangles[j][r];
+ }
+ cell_type_vector[j] = CellType::Triangle;
}
const size_t nb_quadrangles = __quadrangles.extent(0);
for (size_t j=0; j<nb_quadrangles; ++j) {
const size_t jq = j+nb_triangles;
- cell_nodes(jq,0) = __quadrangles[j][0];
- cell_nodes(jq,1) = __quadrangles[j][1];
- cell_nodes(jq,2) = __quadrangles[j][2];
- cell_nodes(jq,3) = __quadrangles[j][3];
- }
- const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
- for (size_t j=0; j<nb_triangles; ++j) {
- cell_nb_nodes[j] = 3;
- }
- for (size_t j=nb_triangles; j<nb_triangles+nb_quadrangles; ++j) {
- cell_nb_nodes[j] = 4;
+ cell_by_node_vector[jq].resize(4);
+ for (int r=0; r<4; ++r) {
+ cell_by_node_vector[jq][r] = __quadrangles[j][r];
+ }
+ cell_type_vector[jq] = CellType::Quadrangle;
}
- std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_nb_nodes, cell_nodes));
+
+ std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector,
+ cell_type_vector));
Connectivity2D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int e=0; e<__edges.extent(0); ++e) {
- const unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
+ const unsigned int edge_number = connectivity.getFaceNumber({__edges[e][0], __edges[e][1]});
const unsigned int& ref = __edges_ref[e];
ref_faces_map[ref].push_back(edge_number);
}
@@ -943,22 +926,21 @@ GmshReader::__proceedData()
m_mesh = std::shared_ptr<IMesh>(new MeshType(p_connectivity, xr));
} else if ((dimension1_mask, elementNumber)>0) {
-
const size_t nb_cells = (dimension1_mask, elementNumber);
- const size_t max_nb_node_per_cell=2;
- const Kokkos::View<unsigned int**> cell_nodes("cell_nodes", nb_cells, max_nb_node_per_cell);
- for (size_t j=0; j<nb_cells; ++j) {
- cell_nodes(j,0) = __edges[j][0];
- cell_nodes(j,1) = __edges[j][1];
- }
+ std::vector<CellType> cell_type_vector(nb_cells);
- const Kokkos::View<unsigned short*> cell_nb_nodes("cell_nb_nodes", nb_cells);
+ std::vector<std::vector<unsigned int>> cell_by_node_vector(nb_cells);
for (size_t j=0; j<nb_cells; ++j) {
- cell_nb_nodes[j] = 2;
+ cell_by_node_vector[j].resize(2);
+ for (int r=0; r<2; ++r) {
+ cell_by_node_vector[j][r] = __edges[j][r];
+ }
+ cell_type_vector[j] = CellType::Line;
}
- std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_nb_nodes, cell_nodes));
+ std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector,
+ cell_type_vector));
Connectivity1D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index f02cfa16f116d2df744cc73bacc073b8964526de..8825a67df3d6cf1593ff54bb6af8b7a96b079b5b 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -4,7 +4,6 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
-#include <Connectivity2D.hpp>
#include <Mesh.hpp>
#include <RefId.hpp>
diff --git a/src/mesh/IConnectivity.hpp b/src/mesh/IConnectivity.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..49432fdea276d770c993960a2ffb48da2f7407b0
--- /dev/null
+++ b/src/mesh/IConnectivity.hpp
@@ -0,0 +1,18 @@
+#ifndef ICONNECTIVITY_HPP
+#define ICONNECTIVITY_HPP
+
+#include <ItemType.hpp>
+#include <ConnectivityMatrix.hpp>
+
+struct IConnectivity
+{
+ virtual const ConnectivityMatrix&
+ getMatrix(const ItemType& item_type_0,
+ const ItemType& item_type_1) const = 0;
+
+ IConnectivity() = default;
+ IConnectivity(const IConnectivity&) = delete;
+ ~IConnectivity() = default;
+};
+
+#endif // ICONNECTIVITY_HPP
diff --git a/src/mesh/ItemType.hpp b/src/mesh/ItemType.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..d37c02fae15120b872c008f939c5d4907debb085
--- /dev/null
+++ b/src/mesh/ItemType.hpp
@@ -0,0 +1,118 @@
+#ifndef ITEM_TYPE_HPP
+#define ITEM_TYPE_HPP
+
+#include <utility>
+#include <limits>
+#include <string>
+
+enum class ItemType
+{
+ node = 0,
+ edge = 1,
+ face = 2,
+ cell = 3
+};
+
+inline constexpr
+std::string_view itemName(const ItemType& item_type)
+{
+ std::string_view name;
+ switch(item_type){
+ case ItemType::node: {
+ name = "node";
+ break;
+ }
+ case ItemType::edge: {
+ name = "edge";
+ break;
+ }
+ case ItemType::face: {
+ name = "face";
+ break;
+ }
+ case ItemType::cell: {
+ name = "cell";
+ break;
+ }
+ }
+ return name;
+}
+
+template <size_t Dimension>
+struct ItemId {};
+
+template <>
+struct ItemId<1>
+{
+ inline static constexpr size_t itemId(const ItemType& item_type) {
+ size_t i = std::numeric_limits<size_t>::max();
+ switch(item_type) {
+ case ItemType::cell: {
+ i=0;
+ break;
+ }
+ case ItemType::edge:
+ case ItemType::face:
+ case ItemType::node: {
+ // in 1d, faces, edges and nodes are the same
+ i=1;
+ break;
+ }
+ }
+ return i;
+ }
+};
+
+template <>
+struct ItemId<2>
+{
+ inline static constexpr size_t itemId(const ItemType& item_type) {
+ size_t i = std::numeric_limits<size_t>::max();
+ switch(item_type) {
+ case ItemType::cell: {
+ i=0;
+ break;
+ }
+ case ItemType::edge:
+ case ItemType::face: {
+ // in 2d, faces and edges are the same
+ i=1;
+ break;
+ }
+ case ItemType::node: {
+ i=2;
+ break;
+ }
+ }
+ return i;
+ }
+};
+
+template <>
+struct ItemId<3>
+{
+ inline static constexpr size_t itemId(const ItemType& item_type) {
+ size_t i = std::numeric_limits<size_t>::max();
+ switch(item_type) {
+ case ItemType::cell: {
+ i=0;
+ break;
+ }
+ case ItemType::edge: {
+ i=1;
+ break;
+ }
+ case ItemType::face: {
+ i=2;
+ break;
+ }
+ case ItemType::node: {
+ i=3;
+ break;
+ }
+ }
+ return i;
+ }
+};
+
+#endif // ITEM_TYPE_HPP
diff --git a/src/mesh/Mesh.hpp b/src/mesh/Mesh.hpp
index 2b05731f54a466635d43996319314968d1da5bb1..3aee2c287d5cdd06c4a754115caf930973c4f2cb 100644
--- a/src/mesh/Mesh.hpp
+++ b/src/mesh/Mesh.hpp
@@ -36,17 +36,20 @@ public:
return *m_connectivity;
}
- const size_t& numberOfNodes() const
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfNodes() const
{
return m_connectivity->numberOfNodes();
}
- const size_t& numberOfFaces() const
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfFaces() const
{
return m_connectivity->numberOfFaces();
}
- const size_t& numberOfCells() const
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfCells() const
{
return m_connectivity->numberOfCells();
}
diff --git a/src/mesh/MeshData.hpp b/src/mesh/MeshData.hpp
index 59411b7f984f318e5af1e35745651208a046abd0..0fe98256bb87d8d02b29899906e475737f4e835d 100644
--- a/src/mesh/MeshData.hpp
+++ b/src/mesh/MeshData.hpp
@@ -4,12 +4,14 @@
#include <Kokkos_Core.hpp>
#include <TinyVector.hpp>
+#include <SubItemValuePerItem.hpp>
+
#include <map>
template <typename M>
class MeshData
{
-public:
+ public:
typedef M MeshType;
static constexpr size_t dimension = MeshType::dimension;
@@ -19,56 +21,62 @@ public:
static constexpr double inv_dimension = 1./dimension;
-private:
+ private:
const MeshType& m_mesh;
- Kokkos::View<Rd**> m_Cjr;
- Kokkos::View<double**> m_ljr;
- Kokkos::View<Rd**> m_njr;
+ NodeValuePerCell<Rd> m_Cjr;
+ NodeValuePerCell<double> m_ljr;
+ NodeValuePerCell<Rd> m_njr;
Kokkos::View<Rd*> m_xj;
Kokkos::View<double*> m_Vj;
KOKKOS_INLINE_FUNCTION
void _updateCenter()
{ // Computes vertices isobarycenter
- if(dimension == 1) {
+ if constexpr (dimension == 1) {
const Kokkos::View<const Rd*> xr = m_mesh.xr();
- const Kokkos::View<const unsigned int**>& cell_nodes
- = m_mesh.connectivity().cellNodes();
+
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- m_xj[j] = 0.5*(xr[cell_nodes(j,0)]+xr[cell_nodes(j,1)]);
- });
+ 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 unsigned int**>& cell_nodes
- = m_mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*>& cell_nb_nodes
- = m_mesh.connectivity().cellNbNodes();
const Kokkos::View<const double*>& inv_cell_nb_nodes
- = m_mesh.connectivity().invCellNbNodes();
+ = m_mesh.connectivity().invCellNbNodes();
+
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- Rd X = zero;
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- X += xr[cell_nodes(j,R)];
- }
- m_xj[j] = inv_cell_nb_nodes[j]*X;
- });
+ 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[cell_nodes(R)];
+ }
+ m_xj[j] = inv_cell_nb_nodes[j]*X;
+ });
}
}
KOKKOS_INLINE_FUNCTION
void _updateVolume()
{
- const Kokkos::View<const unsigned int**>& cell_nodes
- = m_mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_mesh.connectivity().cellNbNodes();
-
const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
double sum_cjr_xr = 0;
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- sum_cjr_xr += (xr[cell_nodes(j,R)], m_Cjr(j,R));
+ 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));
}
m_Vj[j] = inv_dimension * sum_cjr_xr;
});
@@ -78,136 +86,140 @@ private:
void _updateCjr() {
if constexpr (dimension == 1) {
// Cjr/njr/ljr are constant overtime
- }
+ }
else if constexpr (dimension == 2) {
- const Kokkos::View<const unsigned int**>& cell_nodes
- = m_mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_mesh.connectivity().cellNbNodes();
-
const Kokkos::View<const Rd*> xr = m_mesh.xr();
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- int Rp1 = (R+1)%cell_nb_nodes[j];
- int Rm1 = (R+cell_nb_nodes[j]-1)%cell_nb_nodes[j];
- Rd half_xrp_xrm = 0.5*(xr(cell_nodes(j,Rp1))-xr(cell_nodes(j,Rm1)));
+ 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(cell_nodes(Rp1))-xr(cell_nodes(Rm1)));
m_Cjr(j,R) = Rd{-half_xrp_xrm[1], half_xrp_xrm[0]};
}
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const Rd& Cjr = m_Cjr(j,R);
- m_ljr(j,R) = l2Norm(Cjr);
- }
+ const NodeValuePerCell<Rd>& Cjr = m_Cjr;
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
+ m_ljr[j] = l2Norm(Cjr[j]);
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const Rd& Cjr = m_Cjr(j,R);
- const double inv_ljr = 1./m_ljr(j,R);
- m_njr(j,R) = inv_ljr*Cjr;
- }
+ const NodeValuePerCell<double>& ljr = m_ljr;
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& j){
+ m_njr[j] = (1./ljr[j])*Cjr[j];
});
- } else if (dimension ==3) {
- const Kokkos::View<const unsigned int**>& cell_nodes
- = m_mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_mesh.connectivity().cellNbNodes();
+ } else if (dimension ==3) {
const Kokkos::View<const Rd*> xr = m_mesh.xr();
- Kokkos::View<Rd**> Nlr("Nlr", m_mesh.connectivity().numberOfFaces(), m_mesh.connectivity().maxNbNodePerFace());
+ NodeValuePerFace<Rd> Nlr(m_mesh.connectivity());
+ const auto& face_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::face,
+ ItemType::node);
- const Kokkos::View<const unsigned int**>& face_nodes
- = m_mesh.connectivity().faceNodes();
- const Kokkos::View<const unsigned short*>& face_nb_nodes
- = m_mesh.connectivity().faceNbNodes();
Kokkos::parallel_for(m_mesh.numberOfFaces(), KOKKOS_LAMBDA(const int& l) {
- const size_t nb_nodes = face_nb_nodes[l];
+ const auto& face_nodes = face_to_node_matrix.rowConst(l);
+ const size_t nb_nodes = face_nodes.length;
std::vector<Rd> dxr(nb_nodes);
for (size_t r=0; r<nb_nodes; ++r) {
- dxr[r] = xr[face_nodes(l,(r+1)%nb_nodes)] - xr[face_nodes(l,(r+nb_nodes-1)%nb_nodes)];
+ dxr[r] = 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[face_nodes(l,s)]);
+ Nr += crossProduct((two_dxr - dxr[s]), xr[face_nodes(s)]);
}
Nr *= inv_12_nb_nodes;
- Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes(l,r)]);
+ Nr -= (1./6.)*crossProduct(dxr[r], xr[face_nodes(r)]);
Nlr(l,r) = Nr;
}
});
- const Kokkos::View<const unsigned int**> cell_faces
- = m_mesh.connectivity().cellFaces();
- const Kokkos::View<const bool**> cell_faces_is_reversed
- = m_mesh.connectivity().cellFacesIsReversed();
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ m_Cjr[jr] = zero;
+ });
+
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
- const Kokkos::View<const unsigned short*> cell_nb_faces
- = m_mesh.connectivity().cellNbFaces();
+ const auto& cell_to_face_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::face);
+ const auto& cell_face_is_reversed = m_mesh.connectivity().cellFaceIsReversed();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- m_Cjr(j,R) = zero;
- }
- std::map<unsigned int, unsigned short> node_id_to_local;
- for (size_t R=0; R<cell_nb_nodes[j]; ++R) {
- node_id_to_local[cell_nodes(j,R)] = R;
- }
- for (size_t L=0; L<cell_nb_faces[j]; ++L) {
- const size_t l = cell_faces(j, L);
- if (cell_faces_is_reversed(j, L)) {
- for (size_t rl = 0; rl<face_nb_nodes[l]; ++rl) {
- m_Cjr(j, node_id_to_local[face_nodes(l,rl)]) -= Nlr(l, rl);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+
+ const auto& cell_faces = cell_to_face_matrix.rowConst(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);
+
+#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)) {
+ return i_node;
+ break;
+ }
+ }
+ return std::numeric_limits<size_t>::max();
+ };
+
+ 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));
+ m_Cjr(j, R) -= Nlr(l,rl);
}
} else {
- for (size_t rl = 0; rl<face_nb_nodes[l]; ++rl) {
- m_Cjr(j, node_id_to_local[face_nodes(l,rl)]) += Nlr(l, rl);
+ for (size_t rl = 0; rl<face_nodes.length; ++rl) {
+ const size_t R = local_node_number_in_cell(face_nodes(rl));
+ m_Cjr(j, R) += Nlr(l,rl);
}
}
}
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const Rd& Cjr = m_Cjr(j,R);
- m_ljr(j,R) = l2Norm(Cjr);
- }
+ const NodeValuePerCell<Rd>& Cjr = m_Cjr;
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ m_ljr[jr] = l2Norm(Cjr[jr]);
});
- Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const Rd& Cjr = m_Cjr(j,R);
- const double inv_ljr = 1./m_ljr(j,R);
- m_njr(j,R) = inv_ljr*Cjr;
- }
+ const NodeValuePerCell<double>& ljr = m_ljr;
+ Kokkos::parallel_for(m_Cjr.numberOfValues(), KOKKOS_LAMBDA(const int& jr){
+ m_njr[jr] = (1./ljr[jr])*Cjr[jr];
});
}
static_assert((dimension<=3), "only 1d, 2d and 3d are implemented");
}
-public:
+ public:
const MeshType& mesh() const
{
return m_mesh;
}
- const Kokkos::View<const Rd**> Cjr() const
+ const NodeValuePerCell<Rd>& Cjr() const
{
return m_Cjr;
}
- const Kokkos::View<const double**> ljr() const
+ const NodeValuePerCell<double>& ljr() const
{
return m_ljr;
}
- const Kokkos::View<const Rd**> njr() const
+ const NodeValuePerCell<Rd>& njr() const
{
return m_njr;
}
@@ -230,25 +242,25 @@ public:
}
MeshData(const MeshType& mesh)
- : m_mesh(mesh),
- m_Cjr("Cjr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
- m_ljr("ljr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
- m_njr("njr", mesh.numberOfCells(), mesh.connectivity().maxNbNodePerCell()),
- m_xj("xj", mesh.numberOfCells()),
- m_Vj("Vj", mesh.numberOfCells())
+ : m_mesh(mesh),
+ m_Cjr(mesh.connectivity()),
+ m_ljr(mesh.connectivity()),
+ m_njr(mesh.connectivity()),
+ m_xj("xj", mesh.numberOfCells()),
+ m_Vj("Vj", mesh.numberOfCells())
{
if constexpr (dimension==1) {
// in 1d Cjr are computed once for all
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- m_Cjr(j,0)=-1;
- m_Cjr(j,1)= 1;
- });
- // in 1d njr=Cjr
+ m_Cjr(j,0)=-1;
+ m_Cjr(j,1)= 1;
+ });
+ // in 1d njr=Cjr (here performs a shallow copy)
m_njr=m_Cjr;
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- m_ljr(j,0)= 1;
- m_ljr(j,1)= 1;
- });
+ m_ljr(j,0)= 1;
+ m_ljr(j,1)= 1;
+ });
}
this->updateAllData();
}
diff --git a/src/mesh/MeshNodeBoundary.hpp b/src/mesh/MeshNodeBoundary.hpp
index 8361c2f29cab7253f10c76c639016cad54d2fd73..87c537439bf1a27e8443865677cd6b58c132378c 100644
--- a/src/mesh/MeshNodeBoundary.hpp
+++ b/src/mesh/MeshNodeBoundary.hpp
@@ -2,11 +2,14 @@
#define MESH_NODE_BOUNDARY_HPP
#include <Kokkos_Core.hpp>
+#include <Kokkos_Vector.hpp>
#include <TinyVector.hpp>
#include <RefNodeList.hpp>
#include <RefFaceList.hpp>
+#include <ConnectivityMatrix.hpp>
+#include <IConnectivity.hpp>
#include <iostream>
template <size_t dimension>
@@ -28,26 +31,31 @@ class MeshNodeBoundary
const RefFaceList& ref_face_list)
{
static_assert(dimension == MeshType::dimension);
- const Kokkos::View<const unsigned short*> face_nb_cells = mesh.connectivity().faceNbCells();
+ const auto& face_to_cell_matrix
+ = mesh.connectivity().getMatrix(ItemType::face, ItemType::cell);
const Kokkos::View<const unsigned int*>& face_list = ref_face_list.faceList();
Kokkos::parallel_for(face_list.extent(0), KOKKOS_LAMBDA(const int& l){
- if (face_nb_cells[face_list[l]]>1) {
+ const auto& face_cells = face_to_cell_matrix.rowConst(face_list[l]);
+ if (face_cells.length>1) {
std::cerr << "internal faces cannot be used to define mesh boundaries\n";
std::exit(1);
}
});
- const Kokkos::View<const unsigned short*> face_nb_nodes = mesh.connectivity().faceNbNodes();
- const Kokkos::View<const unsigned int**> face_nodes = mesh.connectivity().faceNodes();
-
- std::vector<unsigned int> node_ids;
+ Kokkos::vector<unsigned int> node_ids;
// not enough but should reduce significantly the number of resizing
node_ids.reserve(dimension*face_list.extent(0));
+ const auto& face_to_node_matrix
+ = mesh.connectivity().getMatrix(ItemType::face,
+ ItemType::node);
+
for (size_t l=0; l<face_list.extent(0); ++l) {
const size_t face_number = face_list[l];
- for (size_t r=0; r<face_nb_nodes[face_number]; ++r) {
- node_ids.push_back(face_nodes(face_number,r));
+ const auto& face_nodes = face_to_node_matrix.rowConst(face_number);
+
+ for (size_t r=0; r<face_nodes.length; ++r) {
+ node_ids.push_back(face_nodes(r));
}
}
std::sort(node_ids.begin(), node_ids.end());
@@ -307,17 +315,22 @@ _getOutgoingNormal(const MeshType& mesh)
const R normal = this->_getNormal(mesh);
const Kokkos::View<const R*>& xr = mesh.xr();
- const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells();
- const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes();
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
+ const auto& node_to_cell_matrix
+ = mesh.connectivity().getMatrix(ItemType::node,
+ ItemType::cell);
const size_t r0 = m_node_list[0];
- const size_t j0 = node_cells(r0,0);
+ const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
- for (int r=0; r<cell_nb_nodes(j0); ++r) {
- const double height = (xr(cell_nodes(j0, r))-xr(r0), normal);
+ 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) {
@@ -339,17 +352,22 @@ _getOutgoingNormal(const MeshType& mesh)
const R2 normal = this->_getNormal(mesh);
const Kokkos::View<const R2*>& xr = mesh.xr();
- const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells();
- const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes();
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
+ const auto& node_to_cell_matrix
+ = mesh.connectivity().getMatrix(ItemType::node,
+ ItemType::cell);
const size_t r0 = m_node_list[0];
- const size_t j0 = node_cells(r0,0);
+ const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
- for (int r=0; r<cell_nb_nodes(j0); ++r) {
- const double height = (xr(cell_nodes(j0, r))-xr(r0), normal);
+ 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) {
@@ -371,17 +389,22 @@ _getOutgoingNormal(const MeshType& mesh)
const R3 normal = this->_getNormal(mesh);
const Kokkos::View<const R3*>& xr = mesh.xr();
- const Kokkos::View<const unsigned int**>& node_cells = mesh.connectivity().nodeCells();
- const Kokkos::View<const unsigned int**>& cell_nodes = mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*>& cell_nb_nodes = mesh.connectivity().cellNbNodes();
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
+ const auto& node_to_cell_matrix
+ = mesh.connectivity().getMatrix(ItemType::node,
+ ItemType::cell);
const size_t r0 = m_node_list[0];
- const size_t j0 = node_cells(r0,0);
+ const size_t j0 = node_to_cell_matrix.rowConst(r0)(0);
+ const auto& j0_nodes = cell_to_node_matrix.rowConst(j0);
double max_height = 0;
- for (int r=0; r<cell_nb_nodes(j0); ++r) {
- const double height = (xr(cell_nodes(j0, r))-xr(r0), normal);
+ 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/SubItemValuePerItem.hpp b/src/mesh/SubItemValuePerItem.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..5778fc73989ca547130288d1737bf9b1d767c2ef
--- /dev/null
+++ b/src/mesh/SubItemValuePerItem.hpp
@@ -0,0 +1,249 @@
+#ifndef SUBITEM_VALUE_PER_ITEM_HPP
+#define SUBITEM_VALUE_PER_ITEM_HPP
+
+#include <Kokkos_StaticCrsGraph.hpp>
+#include <ItemType.hpp>
+#include <PastisAssert.hpp>
+
+#include <ConnectivityMatrix.hpp>
+#include <IConnectivity.hpp>
+
+template <typename DataType,
+ ItemType sub_item_type,
+ ItemType item_type,
+ typename Allowed=void>
+class SubItemValuePerItem;
+
+template <typename DataType,
+ ItemType sub_item_type,
+ ItemType item_type>
+class SubItemValuePerItem<DataType,
+ sub_item_type,
+ item_type,
+ std::enable_if_t<sub_item_type != item_type>>
+{
+ public:
+ static const ItemType item_t{item_type};
+ static const ItemType sub_item_t{sub_item_type};
+ private:
+ bool m_is_built{false};
+
+ ConnectivityMatrix::HostRowType m_host_row_map;
+ Kokkos::View<DataType*> m_values;
+
+ // Allows const version to access our data
+ friend SubItemValuePerItem<std::add_const_t<DataType>,
+ sub_item_type,
+ item_type>;
+
+ public:
+ class SubView
+ {
+ private:
+ KOKKOS_RESTRICT DataType* const m_sub_values;
+ const size_t m_size;
+ public:
+ KOKKOS_INLINE_FUNCTION
+ const DataType& operator[](const size_t& i) const
+ {
+ Assert(i<m_size);
+ return m_sub_values[i];
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t& i)
+ {
+ Assert(i<m_size);
+ return m_sub_values[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ const size_t& size() const
+ {
+ return m_size;
+ }
+
+ SubView(const SubView&) = delete;
+
+ KOKKOS_INLINE_FUNCTION
+ SubView(SubView&&) = default;
+
+ KOKKOS_INLINE_FUNCTION
+ SubView(const Kokkos::View<DataType*>& values,
+ const size_t& begin,
+ const size_t& end)
+ : m_sub_values(&(values[begin])),
+ m_size(end-begin)
+ {
+ Assert(begin<=end);
+ Assert(end<=values.extent(0));
+ }
+ };
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ const bool& isBuilt() const
+ {
+ return m_is_built;
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator()(const size_t& j, const size_t& r)
+ {
+ Assert(m_is_built);
+ return m_values[m_host_row_map(j)+r];
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator()(const size_t& j, const size_t& r) const
+ {
+ Assert(m_is_built);
+ return m_values[m_host_row_map(j)+r];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfValues() const
+ {
+ Assert(m_is_built);
+ return m_values.extent(0);
+ }
+
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t& i)
+ {
+ Assert(m_is_built);
+ return m_values[i];
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_FORCEINLINE_FUNCTION
+ DataType& operator[](const size_t & i) const
+ {
+ Assert(m_is_built);
+ return m_values[i];
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfItems() const
+ {
+ Assert(m_is_built);
+ Assert(m_host_row_map.extent(0) != 0);
+ return m_host_row_map.extent(0);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ size_t numberOfSubValues(const size_t& i_cell) const
+ {
+ Assert(m_is_built);
+ return m_host_row_map(i_cell+1)-m_host_row_map(i_cell);
+ }
+
+ KOKKOS_INLINE_FUNCTION
+ SubView itemValues(const size_t& i_cell)
+ {
+ Assert(m_is_built);
+ const auto& cell_begin = m_host_row_map(i_cell);
+ const auto& cell_end = m_host_row_map(i_cell+1);
+ return SubView(m_values, cell_begin, cell_end);
+ }
+
+ // Following Kokkos logic, these classes are view and const view does allow
+ // changes in data
+ KOKKOS_INLINE_FUNCTION
+ SubView itemValues(const size_t& i_cell) const
+ {
+ Assert(m_is_built);
+ const auto& cell_begin = m_host_row_map(i_cell);
+ const auto& cell_end = m_host_row_map(i_cell+1);
+ return SubView(m_values, cell_begin, cell_end);
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ SubItemValuePerItem&
+ operator=(const SubItemValuePerItem<DataType2, sub_item_type, item_type>& sub_item_value_per_item)
+ {
+ // ensures that DataType is the same as source DataType2
+ static_assert(std::is_same<std::remove_const_t<DataType>, std::remove_const_t<DataType2>>(),
+ "Cannot assign SubItemValuePerItem of different type");
+ // ensures that const is not lost through copy
+ static_assert(((std::is_const<DataType2>() and std::is_const<DataType>())
+ or not std::is_const<DataType2>()),
+ "Cannot assign const SubItemValuePerItem to a non const SubItemValuePerItem");
+ m_host_row_map = sub_item_value_per_item.m_host_row_map;
+ m_values = sub_item_value_per_item.m_values;
+
+ m_is_built = sub_item_value_per_item.m_is_built;
+
+ return *this;
+ }
+
+ template <typename DataType2>
+ KOKKOS_INLINE_FUNCTION
+ SubItemValuePerItem(const SubItemValuePerItem<DataType2, sub_item_type, item_type>& sub_item_value_per_item)
+ {
+ this->operator=(sub_item_value_per_item);
+ }
+
+ SubItemValuePerItem() = default;
+
+ SubItemValuePerItem(const IConnectivity& connectivity)
+ : m_is_built{true}
+ {
+ ConnectivityMatrix connectivity_matrix
+ = connectivity.getMatrix(item_type, sub_item_type);
+
+ m_host_row_map = connectivity_matrix.rowsMap();
+ m_values = Kokkos::View<DataType*>("values", connectivity_matrix.numEntries());
+ }
+
+ ~SubItemValuePerItem() = default;
+};
+
+// Item values at nodes
+
+template <typename DataType>
+using NodeValuePerEdge = SubItemValuePerItem<DataType, ItemType::node, ItemType::edge>;
+
+template <typename DataType>
+using NodeValuePerFace = SubItemValuePerItem<DataType, ItemType::node, ItemType::face>;
+
+template <typename DataType>
+using NodeValuePerCell = SubItemValuePerItem<DataType, ItemType::node, ItemType::cell>;
+
+// Item values at edges
+
+template <typename DataType>
+using EdgeValuePerNode = SubItemValuePerItem<DataType, ItemType::edge, ItemType::node>;
+
+template <typename DataType>
+using EdgeValuePerFace = SubItemValuePerItem<DataType, ItemType::edge, ItemType::face>;
+
+template <typename DataType>
+using EdgeValuePerCell = SubItemValuePerItem<DataType, ItemType::edge, ItemType::cell>;
+
+// Item values at faces
+
+template <typename DataType>
+using FaceValuePerNode = SubItemValuePerItem<DataType, ItemType::face, ItemType::node>;
+
+template <typename DataType>
+using FaceValuePerEdge = SubItemValuePerItem<DataType, ItemType::face, ItemType::edge>;
+
+template <typename DataType>
+using FaceValuePerCell = SubItemValuePerItem<DataType, ItemType::face, ItemType::cell>;
+
+// Item values at cells
+
+template <typename DataType>
+using CellValuePerNode = SubItemValuePerItem<DataType, ItemType::cell, ItemType::node>;
+
+template <typename DataType>
+using CellValuePerEdge = SubItemValuePerItem<DataType, ItemType::cell, ItemType::edge>;
+
+template <typename DataType>
+using CellValuePerFace = SubItemValuePerItem<DataType, ItemType::cell, ItemType::face>;
+
+#endif // SUBITEM_VALUE_PER_ITEM_HPP
diff --git a/src/output/VTKWriter.hpp b/src/output/VTKWriter.hpp
index 4e34f666cc0829eef509f4140aaa3aab6d9f9727..a33fef0bfbe2abf2c73b34e7c3b564fb2cd81461 100644
--- a/src/output/VTKWriter.hpp
+++ b/src/output/VTKWriter.hpp
@@ -6,6 +6,7 @@
#include <iomanip>
#include <sstream>
#include <TinyVector.hpp>
+#include <IConnectivity.hpp>
class VTKWriter
{
@@ -67,13 +68,17 @@ class VTKWriter
fout << "</Points>\n";
fout << "<Cells>\n";
- const Kokkos::View<const unsigned int**> cell_nodes = mesh.connectivity().cellNodes();
- const Kokkos::View<const unsigned short*> cell_nb_nodes = mesh.connectivity().cellNbNodes();
fout << "<DataArray type=\"Int32\" Name=\"connectivity\" NumberOfComponents=\"1\" format=\"ascii\">\n";
+
+ const auto& cell_to_node_matrix
+ = mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
+
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- for (unsigned short r=0; r<cell_nb_nodes[j]; ++r) {
- fout << cell_nodes(j,r) << ' ';
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ for (unsigned short r=0; r<cell_nodes.length; ++r) {
+ fout << cell_nodes(r) << ' ';
}
}
fout << '\n';
@@ -83,7 +88,8 @@ class VTKWriter
{
unsigned int offset=0;
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- offset += cell_nb_nodes[j];
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ offset += cell_nodes.length;
fout << offset << ' ';
}
}
@@ -92,7 +98,8 @@ class VTKWriter
fout << "<DataArray type=\"Int8\" Name=\"types\" NumberOfComponents=\"1\" format=\"ascii\">\n";
for (unsigned int j=0; j<mesh.numberOfCells(); ++j) {
- switch (cell_nb_nodes[j]) {
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+ switch (cell_nodes.length) {
case 2: {
fout << "3 ";
break;
diff --git a/src/scheme/AcousticSolver.hpp b/src/scheme/AcousticSolver.hpp
index 3e2d108d7eef1386aa4ec1c2679ecd26c80ec243..0cc6ec0df7a97b91999376873ed7b35a0c5acf6c 100644
--- a/src/scheme/AcousticSolver.hpp
+++ b/src/scheme/AcousticSolver.hpp
@@ -14,6 +14,8 @@
#include <FiniteVolumesEulerUnknowns.hpp>
#include <BoundaryCondition.hpp>
+#include <SubItemValuePerItem.hpp>
+
template<typename MeshData>
class AcousticSolver
{
@@ -30,13 +32,13 @@ class AcousticSolver
typedef TinyVector<dimension> Rd;
typedef TinyMatrix<dimension> Rdd;
-private:
+ private:
struct ReduceMin
{
- private:
+ private:
const Kokkos::View<const double*> x_;
- public:
+ public:
typedef Kokkos::View<const double*>::non_const_value_type value_type;
ReduceMin(const Kokkos::View<const double*>& x) : x_ (x) {}
@@ -47,16 +49,16 @@ private:
void operator() (const size_type i, value_type& update) const
{
if (x_(i) < update) {
- update = x_(i);
+ update = x_(i);
}
}
KOKKOS_INLINE_FUNCTION
void join (volatile value_type& dst,
- const volatile value_type& src) const
+ const volatile value_type& src) const
{
if (src < dst) {
- dst = src;
+ dst = src;
}
}
@@ -67,50 +69,53 @@ private:
}
};
-
KOKKOS_INLINE_FUNCTION
const Kokkos::View<const double*>
computeRhoCj(const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const double*>& cj)
+ const Kokkos::View<const double*>& cj)
{
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- m_rhocj[j] = rhoj[j]*cj[j];
+ m_rhocj[j] = rhoj[j]*cj[j];
});
return m_rhocj;
}
KOKKOS_INLINE_FUNCTION
- const Kokkos::View<const Rdd**>
- computeAjr(const Kokkos::View<const double*>& rhocj,
- const Kokkos::View<const double**>& ljr,
- const Kokkos::View<const Rd**>& njr) {
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_connectivity.cellNbNodes();
-
+ void computeAjr(const Kokkos::View<const double*>& rhocj,
+ const NodeValuePerCell<Rd>& Cjr,
+ const NodeValuePerCell<double>& ljr,
+ const NodeValuePerCell<Rd>& njr)
+ {
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<cell_nb_nodes[j]; ++r) {
- m_Ajr(j,r) = tensorProduct((rhocj(j)*ljr(j,r))*njr(j,r), njr(j,r));
- }
+ const size_t& nb_nodes =m_Ajr.numberOfSubValues(j);
+ const double& rho_c = rhocj(j);
+ for (size_t r=0; r<nb_nodes; ++r) {
+ m_Ajr(j,r) = tensorProduct(rho_c*Cjr(j,r), njr(j,r));
+ }
});
-
- return m_Ajr;
}
KOKKOS_INLINE_FUNCTION
const Kokkos::View<const Rdd*>
- computeAr(const Kokkos::View<const Rdd**>& Ajr) {
- const Kokkos::View<const unsigned int**> node_cells = m_connectivity.nodeCells();
- const Kokkos::View<const unsigned short**> node_cell_local_node = m_connectivity.nodeCellLocalNode();
- const Kokkos::View<const unsigned short*> node_nb_cells = m_connectivity.nodeNbCells();
+ computeAr(const NodeValuePerCell<Rdd>& Ajr) {
+ const auto& node_to_cell_matrix
+ = m_connectivity.getMatrix(ItemType::node,
+ ItemType::cell);
+ const auto& node_local_numbers_in_their_cells
+ = m_connectivity.nodeLocalNumbersInTheirCells();
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- Rdd sum = zero;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- sum += Ajr(J,R);
- }
- m_Ar(r) = sum;
+ Rdd sum = zero;
+ const auto& node_to_cell = node_to_cell_matrix.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 unsigned int J = node_to_cell(j);
+ const unsigned int R = node_local_number_in_its_cells[j];
+ sum += Ajr(J,R);
+ }
+ m_Ar(r) = sum;
});
return m_Ar;
@@ -118,22 +123,27 @@ private:
KOKKOS_INLINE_FUNCTION
const Kokkos::View<const Rd*>
- computeBr(const Kokkos::View<const Rdd**>& Ajr,
- const Kokkos::View<const Rd**>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj) {
- const Kokkos::View<const unsigned int**>& node_cells = m_connectivity.nodeCells();
- const Kokkos::View<const unsigned short**>& node_cell_local_node = m_connectivity.nodeCellLocalNode();
- const Kokkos::View<const unsigned short*>& node_nb_cells = m_connectivity.nodeNbCells();
+ computeBr(const NodeValuePerCell<Rdd>& Ajr,
+ const NodeValuePerCell<Rd>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj) {
+ const auto& node_to_cell_matrix
+ = m_connectivity.getMatrix(ItemType::node,
+ ItemType::cell);
+ const auto& node_local_numbers_in_their_cells
+ = m_connectivity.nodeLocalNumbersInTheirCells();
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- Rd& br = m_br(r);
- br = zero;
- for (int j=0; j<node_nb_cells(r); ++j) {
- const int J = node_cells(r,j);
- const int R = node_cell_local_node(r,j);
- br += Ajr(J,R)*uj(J) + pj(J)*Cjr(J,R);
- }
+ Rd& br = m_br(r);
+ br = zero;
+ const auto& node_to_cell = node_to_cell_matrix.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 unsigned int 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);
+ }
});
return m_br;
@@ -143,171 +153,155 @@ private:
{
for (const auto& handler : m_boundary_condition_list) {
switch (handler.boundaryCondition().type()) {
- case BoundaryCondition::normal_velocity: {
- std::cerr << "normal_velocity NIY\n";
- std::exit(0);
- break;
- }
- case BoundaryCondition::velocity: {
- std::cerr << "velocity NIY\n";
- std::exit(0);
- break;
- }
- case BoundaryCondition::pressure: {
- const PressureBoundaryCondition& pressure_bc
- = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
- const Kokkos::View<const unsigned int**> face_cells
- = m_connectivity.faceCells();
- const Kokkos::View<const unsigned short**> face_cell_local_face
- = m_connectivity.faceCellLocalFace();
- const Kokkos::View<const Rd**> Cjr
- = m_mesh_data.Cjr();
-
- Kokkos::parallel_for(pressure_bc.numberOfFaces(), KOKKOS_LAMBDA(const int& l_number) {
- // quite ugly: node/faces are melt in a bad way... at least works in 1d...
- const int l = pressure_bc.faceList()[l_number];
- Assert(m_connectivity.faceNbCells()[l] == 1);
- const unsigned int j = face_cells(l,0);
- const unsigned int L = face_cell_local_face(l,0);
-
- m_br(l) -= pressure_bc.value()*Cjr(j,L);
- });
- break;
- }
- case BoundaryCondition::symmetry: {
- const SymmetryBoundaryCondition<dimension>& symmetry_bc
- = dynamic_cast<const SymmetryBoundaryCondition<dimension>&>(handler.boundaryCondition());
- const Rd& n = symmetry_bc.outgoingNormal();
-
- const Rdd I = identity;
- const Rdd nxn = tensorProduct(n,n);
- const Rdd P = I-nxn;
-
- Kokkos::parallel_for(symmetry_bc.numberOfNodes(), KOKKOS_LAMBDA(const int& r_number) {
- const int r = symmetry_bc.nodeList()[r_number];
- // assert(m_connectivity.nodeNbCells()[r] == 1);
-
- m_Ar(r) = P*m_Ar(r)*P + nxn;
- m_br(r) = P*m_br(r);
- });
- break;
- }
+ case BoundaryCondition::normal_velocity: {
+ std::cerr << __FILE__ << ':' << __LINE__ << ": normal_velocity BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::velocity: {
+ std::cerr << __FILE__ << ':' << __LINE__ << ": velocity BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::pressure: {
+ // const PressureBoundaryCondition& pressure_bc
+ // = dynamic_cast<const PressureBoundaryCondition&>(handler.boundaryCondition());
+ std::cerr << __FILE__ << ':' << __LINE__ << ": pressure BC NIY\n";
+ std::exit(0);
+ break;
+ }
+ case BoundaryCondition::symmetry: {
+ const SymmetryBoundaryCondition<dimension>& symmetry_bc
+ = dynamic_cast<const SymmetryBoundaryCondition<dimension>&>(handler.boundaryCondition());
+ const Rd& n = symmetry_bc.outgoingNormal();
+
+ const Rdd I = identity;
+ const Rdd nxn = tensorProduct(n,n);
+ const Rdd P = I-nxn;
+
+ Kokkos::parallel_for(symmetry_bc.numberOfNodes(), KOKKOS_LAMBDA(const int& r_number) {
+ const int r = symmetry_bc.nodeList()[r_number];
+
+ m_Ar(r) = P*m_Ar(r)*P + nxn;
+ m_br(r) = P*m_br(r);
+ });
+ break;
+ }
}
}
}
Kokkos::View<Rd*>
computeUr(const Kokkos::View<const Rdd*>& Ar,
- const Kokkos::View<const Rd*>& br) {
+ const Kokkos::View<const Rd*>& br)
+ {
inverse(Ar, m_inv_Ar);
const Kokkos::View<const Rdd*> invAr = m_inv_Ar;
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
- m_ur[r]=invAr(r)*br(r);
+ m_ur[r]=invAr(r)*br(r);
});
return m_ur;
}
- Kokkos::View<Rd**>
- computeFjr(const Kokkos::View<const Rdd**>& Ajr,
- const Kokkos::View<const Rd*>& ur,
- const Kokkos::View<const Rd**>& Cjr,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj) {
- const Kokkos::View<const unsigned int**>& cell_nodes = m_connectivity.cellNodes();
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_connectivity.cellNbNodes();
+ void
+ computeFjr(const NodeValuePerCell<Rdd>& Ajr,
+ const Kokkos::View<const Rd*>& ur,
+ const NodeValuePerCell<Rd>& Cjr,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj)
+ {
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- for (int r=0; r<cell_nb_nodes[j]; ++r) {
- m_Fjr(j,r) = Ajr(j,r)*(uj(j)-ur(cell_nodes(j,r)))+pj(j)*Cjr(j,r);
- }
- });
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
- return m_Fjr;
+ 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);
+ }
+ });
}
void inverse(const Kokkos::View<const Rdd*>& A,
- Kokkos::View<Rdd*>& inv_A) const {
+ Kokkos::View<Rdd*>& inv_A) const
+ {
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r) {
inv_A(r) = ::inverse(A(r));
});
}
- void inverse(const Kokkos::View<const double*>& x,
- Kokkos::View<double*>& inv_x) const {
- Kokkos::parallel_for(x.size(), KOKKOS_LAMBDA(const int& r) {
- inv_x(r) = 1./x(r);
- });
- }
-
KOKKOS_INLINE_FUNCTION
void computeExplicitFluxes(const Kokkos::View<const Rd*>& xr,
- const Kokkos::View<const Rd*>& xj,
- const Kokkos::View<const double*>& rhoj,
- const Kokkos::View<const Rd*>& uj,
- const Kokkos::View<const double*>& pj,
- const Kokkos::View<const double*>& cj,
- const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const Rd**>& Cjr,
- const Kokkos::View<const double**>& ljr,
- const Kokkos::View<const Rd**>& njr) {
+ const Kokkos::View<const Rd*>& xj,
+ const Kokkos::View<const double*>& rhoj,
+ const Kokkos::View<const Rd*>& uj,
+ const Kokkos::View<const double*>& pj,
+ const Kokkos::View<const double*>& cj,
+ const Kokkos::View<const double*>& Vj,
+ const NodeValuePerCell<Rd>& Cjr,
+ const NodeValuePerCell<double>& ljr,
+ const NodeValuePerCell<Rd>& njr) {
const Kokkos::View<const double*> rhocj = computeRhoCj(rhoj, cj);
- const Kokkos::View<const Rdd**> Ajr = computeAjr(rhocj, ljr, njr);
+ computeAjr(rhocj, Cjr, ljr, njr);
- const Kokkos::View<const Rdd*> Ar = computeAr(Ajr);
- const Kokkos::View<const Rd*> br = computeBr(Ajr, Cjr, uj, pj);
+ const Kokkos::View<const Rdd*> Ar = computeAr(m_Ajr);
+ const Kokkos::View<const Rd*> br = computeBr(m_Ajr, Cjr, uj, pj);
this->applyBoundaryConditions();
Kokkos::View<Rd*> ur = m_ur;
- Kokkos::View<Rd**> Fjr = m_Fjr;
ur = computeUr(Ar, br);
- Fjr = computeFjr(Ajr, ur, Cjr, uj, pj);
+ computeFjr(m_Ajr, ur, Cjr, uj, pj);
}
Kokkos::View<Rd*> m_br;
- Kokkos::View<Rdd**> m_Ajr;
+ NodeValuePerCell<Rdd> m_Ajr;
Kokkos::View<Rdd*> m_Ar;
Kokkos::View<Rdd*> m_inv_Ar;
- Kokkos::View<Rd**> m_Fjr;
+ NodeValuePerCell<Rd> m_Fjr;
Kokkos::View<Rd*> m_ur;
Kokkos::View<double*> m_rhocj;
Kokkos::View<double*> m_Vj_over_cj;
-public:
+ public:
AcousticSolver(MeshData& mesh_data,
- UnknownsType& unknowns,
- const std::vector<BoundaryConditionHandler>& bc_list)
- : m_mesh_data(mesh_data),
- m_mesh(mesh_data.mesh()),
- m_connectivity(m_mesh.connectivity()),
- m_boundary_condition_list(bc_list),
- m_br("br", m_mesh.numberOfNodes()),
- m_Ajr("Ajr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
- m_Ar("Ar", m_mesh.numberOfNodes()),
- m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
- m_Fjr("Fjr", m_mesh.numberOfCells(), m_connectivity.maxNbNodePerCell()),
- m_ur("ur", m_mesh.numberOfNodes()),
- m_rhocj("rho_c", m_mesh.numberOfCells()),
- m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
+ UnknownsType& unknowns,
+ const std::vector<BoundaryConditionHandler>& bc_list)
+ : m_mesh_data(mesh_data),
+ m_mesh(mesh_data.mesh()),
+ m_connectivity(m_mesh.connectivity()),
+ m_boundary_condition_list(bc_list),
+ m_br("br", m_mesh.numberOfNodes()),
+ m_Ajr(m_connectivity),
+ m_Ar("Ar", m_mesh.numberOfNodes()),
+ m_inv_Ar("inv_Ar", m_mesh.numberOfNodes()),
+ m_Fjr(m_connectivity),
+ m_ur("ur", m_mesh.numberOfNodes()),
+ m_rhocj("rho_c", m_mesh.numberOfCells()),
+ m_Vj_over_cj("Vj_over_cj", m_mesh.numberOfCells())
{
;
}
KOKKOS_INLINE_FUNCTION
double acoustic_dt(const Kokkos::View<const double*>& Vj,
- const Kokkos::View<const double*>& cj) const {
- const Kokkos::View<const double**> ljr = m_mesh_data.ljr();
- const Kokkos::View<const unsigned short*>& cell_nb_nodes
- = m_connectivity.cellNbNodes();
+ const Kokkos::View<const double*>& cj) const
+ {
+ const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- double S = 0;
- for (int r=0; r<cell_nb_nodes(j); ++r) {
- S += ljr(j,r);
- }
- m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+
+ double S = 0;
+ for (size_t r=0; r<cell_nodes.length; ++r) {
+ S += ljr(j,r);
+ }
+ m_Vj_over_cj[j] = 2*Vj[j]/(S*cj[j]);
});
double dt = std::numeric_limits<double>::max();
@@ -318,7 +312,7 @@ public:
void computeNextStep(const double& t, const double& dt,
- UnknownsType& unknowns)
+ UnknownsType& unknowns)
{
Kokkos::View<double*> rhoj = unknowns.rhoj();
Kokkos::View<Rd*> uj = unknowns.uj();
@@ -331,46 +325,47 @@ public:
const Kokkos::View<const Rd*> xj = m_mesh_data.xj();
const Kokkos::View<const double*> Vj = m_mesh_data.Vj();
- const Kokkos::View<const Rd**> Cjr = m_mesh_data.Cjr();
- const Kokkos::View<const double**> ljr = m_mesh_data.ljr();
- const Kokkos::View<const Rd**> njr = m_mesh_data.njr();
+ const NodeValuePerCell<Rd>& Cjr = m_mesh_data.Cjr();
+ const NodeValuePerCell<double>& ljr = m_mesh_data.ljr();
+ const NodeValuePerCell<Rd>& njr = m_mesh_data.njr();
Kokkos::View<Rd*> xr = m_mesh.xr();
computeExplicitFluxes(xr, xj, rhoj, uj, pj, cj, Vj, Cjr, ljr, njr);
- const Kokkos::View<const Rd**> Fjr = m_Fjr;
+ const NodeValuePerCell<Rd>& Fjr = m_Fjr;
const Kokkos::View<const Rd*> ur = m_ur;
- const Kokkos::View<const unsigned short*> cell_nb_nodes
- = m_connectivity.cellNbNodes();
- const Kokkos::View<const unsigned int**>& cell_nodes
- = m_connectivity.cellNodes();
+ const auto& cell_to_node_matrix
+ = m_mesh.connectivity().getMatrix(ItemType::cell,
+ ItemType::node);
const Kokkos::View<const double*> inv_mj = unknowns.invMj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- Rd momentum_fluxes = zero;
- double energy_fluxes = 0;
- for (int R=0; R<cell_nb_nodes[j]; ++R) {
- const int r=cell_nodes(j,R);
- momentum_fluxes += Fjr(j,R);
- energy_fluxes += (Fjr(j,R), ur[r]);
- }
- uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
- Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
+ const auto& cell_nodes = cell_to_node_matrix.rowConst(j);
+
+ Rd momentum_fluxes = zero;
+ double energy_fluxes = 0;
+ for (size_t R=0; R<cell_nodes.length; ++R) {
+ const unsigned int r=cell_nodes(R);
+ momentum_fluxes += Fjr(j,R);
+ energy_fluxes += (Fjr(j,R), ur[r]);
+ }
+ uj[j] -= (dt*inv_mj[j]) * momentum_fluxes;
+ Ej[j] -= (dt*inv_mj[j]) * energy_fluxes;
});
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j) {
- ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
+ ej[j] = Ej[j] - 0.5 * (uj[j],uj[j]);
});
Kokkos::parallel_for(m_mesh.numberOfNodes(), KOKKOS_LAMBDA(const int& r){
- xr[r] += dt*ur[r];
+ xr[r] += dt*ur[r];
});
m_mesh_data.updateAllData();
const Kokkos::View<const double*> mj = unknowns.mj();
Kokkos::parallel_for(m_mesh.numberOfCells(), KOKKOS_LAMBDA(const int& j){
- rhoj[j] = mj[j]/Vj[j];
+ rhoj[j] = mj[j]/Vj[j];
});
}
};
diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt
index 2c1ac56dc2e5ba1886c5efc79709ea6f01deecd4..3b4823f1eaa5d56c7bb2a7e7201d623e8308d60a 100644
--- a/tests/CMakeLists.txt
+++ b/tests/CMakeLists.txt
@@ -4,6 +4,7 @@ set(EXECUTABLE_OUTPUT_PATH ${PASTIS_BINARY_DIR})
add_executable (unit_tests
test_main.cpp
+ test_ItemType.cpp
test_PastisAssert.cpp
test_RevisionInfo.cpp
test_TinyMatrix.cpp
diff --git a/tests/test_ItemType.cpp b/tests/test_ItemType.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..7c39963613dd3bae8d60e5fb4c9a4b0915ce6343
--- /dev/null
+++ b/tests/test_ItemType.cpp
@@ -0,0 +1,48 @@
+#include <catch.hpp>
+
+#include <ItemType.hpp>
+
+TEST_CASE("ItemType", "[connectivity]") {
+
+ ItemType node_type = ItemType::node;
+ ItemType edge_type = ItemType::edge;
+ ItemType face_type = ItemType::face;
+ ItemType cell_type = ItemType::cell;
+
+ SECTION("checking for item type differences") {
+ REQUIRE(((node_type != edge_type) and
+ (node_type != face_type) and
+ (node_type != cell_type) and
+ (edge_type != face_type) and
+ (edge_type != cell_type) and
+ (face_type != cell_type)));
+ }
+
+ SECTION("checking for item type names") {
+ REQUIRE(itemName(node_type)=="node");
+ REQUIRE(itemName(edge_type)=="edge");
+ REQUIRE(itemName(face_type)=="face");
+ REQUIRE(itemName(cell_type)=="cell");
+ }
+
+ SECTION("checking for item ids in 1d") {
+ REQUIRE(ItemId<1>::itemId(cell_type)==0);
+ REQUIRE(ItemId<1>::itemId(edge_type)==1);
+ REQUIRE(ItemId<1>::itemId(face_type)==1);
+ REQUIRE(ItemId<1>::itemId(node_type)==1);
+ }
+
+ SECTION("checking for item ids in 2d") {
+ REQUIRE(ItemId<2>::itemId(cell_type)==0);
+ REQUIRE(ItemId<2>::itemId(edge_type)==1);
+ REQUIRE(ItemId<2>::itemId(face_type)==1);
+ REQUIRE(ItemId<2>::itemId(node_type)==2);
+ }
+
+ SECTION("checking for item ids in 3d") {
+ REQUIRE(ItemId<3>::itemId(cell_type)==0);
+ REQUIRE(ItemId<3>::itemId(edge_type)==1);
+ REQUIRE(ItemId<3>::itemId(face_type)==2);
+ REQUIRE(ItemId<3>::itemId(node_type)==3);
+ }
+}
diff --git a/tests/test_PastisAssert.cpp b/tests/test_PastisAssert.cpp
index 8d07321a113be496a220cbfe62934248c3c8907b..d55ad8d826d75ac10e81b5b110495e74d131e921 100644
--- a/tests/test_PastisAssert.cpp
+++ b/tests/test_PastisAssert.cpp
@@ -3,7 +3,7 @@
#include <PastisAssert.hpp>
#include <string>
-TEST_CASE("PastisAssert", "[assert]") {
+TEST_CASE("PastisAssert", "[utils]") {
SECTION("checking for assert error") {
const std::string filename = "filename";
const int line = 10;