diff --git a/src/mesh/Connectivity2D.hpp b/src/mesh/Connectivity2D.hpp
index 66603aaa86cc0546c5657eb11f2793d35868bb5a..fe7eb44d09612d2ca1ec0613b38ae160b08c9ba2 100644
--- a/src/mesh/Connectivity2D.hpp
+++ b/src/mesh/Connectivity2D.hpp
@@ -22,7 +22,7 @@ private:
const Kokkos::View<const unsigned int**> m_cell_nodes;
Kokkos::View<double*> m_inv_cell_nb_nodes;
- Kokkos::View<unsigned short*> m_cell_nb_faces;
+ Kokkos::View<const unsigned short*> m_cell_nb_faces;
Kokkos::View<unsigned int**> m_cell_faces;
Kokkos::View<unsigned short*> m_node_nb_cells;
@@ -33,8 +33,222 @@ private:
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;
+ void _computeNodeCellConnectivities()
+ {
+ // Computes inefficiently node->cells connectivity [Version 0]
+ std::multimap<unsigned int, unsigned int> node_cells_map;
+ using namespace Kokkos::Experimental;
+ Kokkos::parallel_reduce(m_number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
+ const size_t n = m_cell_nb_nodes[j];
+ if (n > nb_max) nb_max = n;
+ }, Max<size_t>(m_max_nb_node_per_cell));
+
+ for (unsigned int j=0; j<m_number_of_cells; ++j) {
+ for (unsigned int r=0; r<m_cell_nb_nodes[j]; ++r) {
+ node_cells_map.insert(std::make_pair(m_cell_nodes(j,r),j));
+ }
+ }
+
+ std::vector<unsigned int> node_ids;
+ node_ids.reserve(node_cells_map.size());
+ for (const auto& node_cell: node_cells_map) {
+ node_ids.push_back(node_cell.first);
+ }
+ auto last_unique = std::unique(node_ids.begin(), node_ids.end());
+ node_ids.resize(std::distance(node_ids.begin(), last_unique));
+
+ m_number_of_nodes = node_ids.size();
+
+ if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) {
+ std::cerr << "sparse node numerotation NIY\n";
+ for (size_t i=0; i<node_ids.size(); ++i) {
+ std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
+ }
+ std::exit(0);
+ }
+
+ std::vector<std::vector<unsigned int>> node_cells_vector(node_ids.size());
+ for (const auto& node_cell: node_cells_map) {
+ node_cells_vector[node_cell.first].push_back(node_cell.second);
+ }
+
+ Kokkos::View<unsigned short*> node_nb_cells("node_nb_cells", node_ids.size());
+ size_t max_node_cells = 0;
+ for (size_t i=0; i<node_cells_vector.size(); ++i) {
+ const auto& cells_vector = node_cells_vector[i];
+ const size_t nb_cells = cells_vector.size();
+ node_nb_cells[i] = nb_cells;
+ if (nb_cells > max_node_cells) {
+ max_node_cells = nb_cells;
+ }
+ }
+ m_node_nb_cells = node_nb_cells;
+ Kokkos::View<unsigned int**> node_cells("node_cells", node_ids.size(), max_node_cells);
+ for (size_t i=0; i<node_cells_vector.size(); ++i) {
+ const auto& cells_vector = node_cells_vector[i];
+ for (size_t j=0; j<cells_vector.size(); ++j) {
+ node_cells(i,j) = cells_vector[j];
+ }
+ }
+ m_node_cells = node_cells;
+
+ Kokkos::View<unsigned short**> node_cell_local_node("node_cell_local_node",
+ node_ids.size(), max_node_cells);
+ Kokkos::parallel_for(m_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<m_cell_nb_nodes[j]; ++R) {
+ if (m_cell_nodes(j,R) == r) {
+ node_cell_local_node(r,J)=R;
+ break;
+ }
+ }
+ }
+ });
+
+ m_node_cell_local_node = node_cell_local_node;
+ }
+
+ 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();
+ std::cout << "number of faces: " << m_number_of_faces << '\n';
+ 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*[2]> 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:
const size_t& numberOfNodes() const
{
@@ -111,8 +325,14 @@ public:
return m_face_cell_local_face;
}
+ unsigned int getFaceNumber(const unsigned int node0_id,
+ const unsigned int node1_id) const
+ {
+ 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)),
@@ -128,81 +348,9 @@ public:
}
assert(m_number_of_cells>0);
- // Computes inefficiently node->cells connectivity [Version 0]
-
- std::multimap<unsigned int, unsigned int> node_cells_map;
- using namespace Kokkos::Experimental;
- Kokkos::parallel_reduce(m_number_of_cells, KOKKOS_LAMBDA(const int& j, size_t& nb_max) {
- const size_t n = m_cell_nb_nodes[j];
- if (n > nb_max) nb_max = n;
- }, Max<size_t>(m_max_nb_node_per_cell));
-
- for (unsigned int j=0; j<m_number_of_cells; ++j) {
- for (unsigned int r=0; r<m_cell_nb_nodes[j]; ++r) {
- node_cells_map.insert(std::make_pair(cell_nodes(j,r),j));
- }
- }
-
- std::vector<unsigned int> node_ids;
- node_ids.reserve(node_cells_map.size());
- for (const auto& node_cell: node_cells_map) {
- node_ids.push_back(node_cell.first);
- }
- auto last_unique = std::unique(node_ids.begin(), node_ids.end());
- node_ids.resize(std::distance(node_ids.begin(), last_unique));
-
- m_number_of_nodes = node_ids.size();
-
- std::cout << "node_ids.size()=" << node_ids.size() << '\n';
- if ((node_ids[0] != 0) or (node_ids[node_ids.size()-1] != node_ids.size()-1)) {
- std::cerr << "sparse node numerotation NIY\n";
- for (size_t i=0; i<node_ids.size(); ++i) {
- std::cout << "node_ids[" << i << "] = " << node_ids[i] << '\n';
- }
- std::exit(0);
- }
-
- std::vector<std::vector<unsigned int>> node_cells_vector(node_ids.size());
- for (const auto& node_cell: node_cells_map) {
- node_cells_vector[node_cell.first].push_back(node_cell.second);
- }
-
-
- Kokkos::View<unsigned short*> node_nb_cells("node_nb_cells", node_ids.size());
- size_t max_node_cells = 0;
- for (size_t i=0; i<node_cells_vector.size(); ++i) {
- const auto& cells_vector = node_cells_vector[i];
- const size_t nb_cells = cells_vector.size();
- node_nb_cells[i] = nb_cells;
- if (nb_cells > max_node_cells) {
- max_node_cells = nb_cells;
- }
- }
- m_node_nb_cells = node_nb_cells;
- Kokkos::View<unsigned int**> node_cells("node_cells", node_ids.size(), max_node_cells);
- for (size_t i=0; i<node_cells_vector.size(); ++i) {
- const auto& cells_vector = node_cells_vector[i];
- for (size_t j=0; j<cells_vector.size(); ++j) {
- node_cells(i,j) = cells_vector[j];
- }
- }
- m_node_cells = node_cells;
-
- Kokkos::View<unsigned short**> node_cell_local_node("node_cell_local_node",
- node_ids.size(), max_node_cells);
- Kokkos::parallel_for(m_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) {
- node_cell_local_node(r,J)=R;
- break;
- }
- }
- }
- });
-
- m_node_cell_local_node = node_cell_local_node;
+ this->_computeNodeCellConnectivities();
+ this->_computeFaceCellConnectivities();
+ //this->_computeNodeFaceConnectivities();
}
~Connectivity2D()
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 6a6268f24e447fe8144b09d7d5296a7358fc9f84..6577a60bd9587644b47324bf0dca404cfe5dfb91 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -799,6 +799,15 @@ GmshReader::__proceedData()
}
m_connectivity = new Connectivity2D(cell_nb_nodes, cell_nodes);
Connectivity2D& connectivity = *m_connectivity;
+
+ for (unsigned int e=0; e<__edges.extent(0); ++e) {
+ unsigned int edge_number = connectivity.getFaceNumber(__edges[e][0], __edges[e][1]);
+ }
+ // __edges[edgeNumber]
+ // = Edge(a, b);
+ // __edges_ref[edgeNumber] = __references[i];
+
+
typedef Mesh<Connectivity2D> MeshType;
typedef TinyVector<2, double> Rd;
@@ -809,6 +818,7 @@ GmshReader::__proceedData()
}
m_mesh = new MeshType(connectivity, xr);
MeshType& mesh = *m_mesh;
+
std::ofstream gnuplot("mesh.gnu");
for (size_t j=0; j<mesh.numberOfCells(); ++j) {