From a91b700d9fbbbb60ef156c7e62be6e833bc4d93d Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 23 Jan 2019 15:03:01 +0100
Subject: [PATCH] Reorient faces following the increasing order of their node
numbers
Previously they were oriented according to the increasing order of the ids of
their nodes. Since NodeId's are local identifiers in a given partition, it
prohibited reproducible parallel calculations independently of the number of
processors.
---
src/mesh/Connectivity.cpp | 32 ++---
src/mesh/GmshReader.cpp | 243 +++++++++++++++++++-------------------
src/mesh/GmshReader.hpp | 5 +-
3 files changed, 137 insertions(+), 143 deletions(-)
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 7c7f7c726..4ec8145a8 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -72,29 +72,21 @@ Connectivity(const ConnectivityDescriptor& descriptor)
= m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
face_to_node_matrix = descriptor.face_to_node_vector;
- if constexpr (Dimension==2) {
- auto& face_to_cell_matrix
- = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::cell)];
- face_to_cell_matrix = descriptor.face_to_cell_vector;
-
- } else if constexpr (Dimension==3) {
- auto& cell_to_face_matrix
- = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
- cell_to_face_matrix = descriptor.cell_to_face_vector;
-
- {
- FaceValuePerCell<bool> cell_face_is_reversed(*this);
- for (CellId j=0; j<descriptor.cell_face_is_reversed_vector.size(); ++j) {
- const auto& face_cells_vector = descriptor.cell_face_is_reversed_vector[j];
- for (unsigned short lj=0; lj<face_cells_vector.size(); ++lj) {
- cell_face_is_reversed(j,lj) = face_cells_vector[lj];
- }
- }
+ auto& cell_to_face_matrix
+ = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)];
+ cell_to_face_matrix = descriptor.cell_to_face_vector;
- m_cell_face_is_reversed = cell_face_is_reversed;
+ {
+ FaceValuePerCell<bool> cell_face_is_reversed(*this);
+ for (CellId j=0; j<descriptor.cell_face_is_reversed_vector.size(); ++j) {
+ const auto& face_cells_vector = descriptor.cell_face_is_reversed_vector[j];
+ for (unsigned short lj=0; lj<face_cells_vector.size(); ++lj) {
+ cell_face_is_reversed(j,lj) = face_cells_vector[lj];
+ }
}
- }
+ m_cell_face_is_reversed = cell_face_is_reversed;
+ }
{
FaceValue<int> face_number(*this);
face_number = convert_to_array(descriptor.face_number_vector);
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index fda25db04..68e4dd3c4 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -132,17 +132,6 @@ ErrorHandler(const std::string& filename,
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>
{
@@ -158,13 +147,24 @@ class ConnectivityFace<2>
}
};
+ private:
+ const std::vector<int>& m_node_number_vector;
+
unsigned int m_node0_id;
unsigned int m_node1_id;
- friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
+ bool m_reversed;
+
+ public:
+
+ std::vector<unsigned int> nodeIdList() const
{
- os << f.m_node0_id << ' ' << f.m_node1_id << ' ';
- return os;
+ return {m_node0_id, m_node1_id};
+ }
+
+ bool reversed() const
+ {
+ return m_reversed;
}
PASTIS_INLINE
@@ -177,9 +177,9 @@ class ConnectivityFace<2>
PASTIS_INLINE
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)));
+ return ((m_node_number_vector[m_node0_id] < m_node_number_vector[f.m_node0_id]) or
+ ((m_node_number_vector[m_node0_id] == m_node_number_vector[f.m_node0_id]) and
+ (m_node_number_vector[m_node1_id]<m_node_number_vector[f.m_node1_id])));
}
PASTIS_INLINE
@@ -189,13 +189,21 @@ class ConnectivityFace<2>
ConnectivityFace& operator=(ConnectivityFace&&) = default;
PASTIS_INLINE
- ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
+ ConnectivityFace(const std::vector<unsigned int>& node_id_list,
+ const std::vector<int>& node_number_vector)
+ : m_node_number_vector(node_number_vector)
{
- 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;
+ Assert(node_id_list.size()==2);
+
+ if (m_node_number_vector[node_id_list[0]] < m_node_number_vector[node_id_list[1]]) {
+ m_node0_id = node_id_list[0];
+ m_node1_id = node_id_list[1];
+ m_reversed = false;
+ } else {
+ m_node0_id = node_id_list[1];
+ m_node1_id = node_id_list[0];
+ m_reversed = true;
+ }
}
PASTIS_INLINE
@@ -225,28 +233,10 @@ class ConnectivityFace<3>
}
};
+ private:
bool m_reversed;
std::vector<NodeId::base_type> m_node_id_list;
-
- friend std::ostream& operator<<(std::ostream& os, const ConnectivityFace& f)
- {
- for (auto id : f.m_node_id_list) {
- os << id << ' ';
- }
- return os;
- }
-
- PASTIS_INLINE
- const bool& reversed() const
- {
- return m_reversed;
- }
-
- PASTIS_INLINE
- const std::vector<unsigned int>& nodeIdList() const
- {
- return m_node_id_list;
- }
+ const std::vector<int>& m_node_number_vector;
PASTIS_INLINE
std::vector<unsigned int> _sort(const std::vector<unsigned int>& node_list)
@@ -269,10 +259,25 @@ class ConnectivityFace<3>
return rotated_node_list;
}
+ public:
+ PASTIS_INLINE
+ const bool& reversed() const
+ {
+ return m_reversed;
+ }
+
+ PASTIS_INLINE
+ const std::vector<unsigned int>& nodeIdList() const
+ {
+ return m_node_id_list;
+ }
+
PASTIS_INLINE
- ConnectivityFace(const std::vector<unsigned int>& given_node_id_list)
+ ConnectivityFace(const std::vector<unsigned int>& given_node_id_list,
+ const std::vector<int>& node_number_vector)
: m_reversed(false),
- m_node_id_list(_sort(given_node_id_list))
+ m_node_id_list(_sort(given_node_id_list)),
+ m_node_number_vector(node_number_vector)
{
;
}
@@ -1168,98 +1173,89 @@ __readPhysicalNames2_2()
}
}
+template <size_t Dimension>
void
-GmshReader::__compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
- // const auto& cell_to_node_matrix
- // = this->_getMatrix(ItemType::cell, ItemType::node);
- using Face = ConnectivityFace<2>;
+ static_assert((Dimension==2) or (Dimension==3),
+ "Invalid dimension to compute cell-face connectivities");
+ using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
+ using Face = ConnectivityFace<Dimension>;
- // In 2D faces are simply define
- using CellFaceId = std::pair<CellId, unsigned short>;
- std::map<Face, std::vector<CellFaceId>> face_cells_map;
+ const auto& node_number_vector = descriptor.node_number_vector;
+ Array<unsigned short> cell_nb_faces(descriptor.cell_by_node_vector.size());
+ std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
const auto& cell_nodes = descriptor.cell_by_node_vector[j];
- for (unsigned short r=0; r<cell_nodes.size(); ++r) {
- NodeId node0_id = cell_nodes[r];
- NodeId node1_id = cell_nodes[(r+1)%cell_nodes.size()];
- 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));
- }
- }
- {
- descriptor.face_to_node_vector.resize(face_cells_map.size());
- int l=0;
- for (const auto& face_info : face_cells_map) {
- const Face& face = face_info.first;
- descriptor.face_to_node_vector[l] = {face.m_node0_id, face.m_node1_id};
- ++l;
- }
- }
-
- {
- descriptor.face_to_cell_vector.resize(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) {
- descriptor.face_to_cell_vector[l].push_back(cell_info.first);
- }
- ++l;
- }
- }
+ if constexpr (Dimension==2) {
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Triangle: {
+ cell_nb_faces[j] = 3;
+ // face 0
+ Face f0({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
- {
- descriptor.face_number_vector.resize(face_cells_map.size());
- for (size_t l=0; l<face_cells_map.size(); ++l) {
- descriptor.face_number_vector[l] = l;
- }
- perr() << __FILE__ << ':' << __LINE__ << ": " << rang::fgB::red << "must use provided faces number" << rang::fg::reset << '\n';
- }
-}
+ // face 1
+ Face f1({cell_nodes[2], cell_nodes[0]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
-void
-GmshReader::__compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
-{
- using CellFaceInfo = std::tuple<CellId, unsigned short, bool>;
- using Face = ConnectivityFace<3>;
+ // face 2
+ Face f2({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+ break;
+ }
+ case CellType::Quadrangle: {
+ cell_nb_faces[j] = 4;
+ // face 0
+ Face f0({cell_nodes[0], cell_nodes[1]}, node_number_vector);
+ face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
- // const auto& cell_to_node_matrix
- // = this->_getMatrix(ItemType::cell, ItemType::node);
+ // face 1
+ Face f1({cell_nodes[1], cell_nodes[2]}, node_number_vector);
+ face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
- Array<unsigned short> cell_nb_faces(descriptor.cell_by_node_vector.size());
- std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
- for (CellId j=0; j<descriptor.cell_by_node_vector.size(); ++j) {
- const auto& cell_nodes = descriptor.cell_by_node_vector[j];
+ // face 2
+ Face f2({cell_nodes[2], cell_nodes[3]}, node_number_vector);
+ face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
- switch (descriptor.cell_type_vector[j]) {
+ // face 3
+ Face f3({cell_nodes[3], cell_nodes[0]}, node_number_vector);
+ face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ break;
+ }
+ default: {
+ perr() << name(descriptor.cell_type_vector[j])
+ << ": unexpected cell type in dimension 2!\n";
+ std::terminate();
+ }
+ }
+ } else if constexpr (Dimension==3) {
+ switch (descriptor.cell_type_vector[j]) {
case CellType::Tetrahedron: {
cell_nb_faces[j] = 4;
// face 0
Face f0({cell_nodes[1],
cell_nodes[2],
- cell_nodes[3]});
+ cell_nodes[3]}, node_number_vector);
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]});
+ cell_nodes[2]}, node_number_vector);
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]});
+ cell_nodes[3]}, node_number_vector);
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]});
+ cell_nodes[1]}, node_number_vector);
face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
break;
}
@@ -1268,50 +1264,52 @@ GmshReader::__compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor&
Face f0({cell_nodes[3],
cell_nodes[2],
cell_nodes[1],
- cell_nodes[0]});
+ cell_nodes[0]}, node_number_vector);
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]});
+ cell_nodes[7]}, node_number_vector);
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]});
+ cell_nodes[3]}, node_number_vector);
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]});
+ cell_nodes[5]}, node_number_vector);
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]});
+ cell_nodes[4]}, node_number_vector);
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]});
+ cell_nodes[2]}, node_number_vector);
face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
cell_nb_faces[j] = 6;
break;
}
default: {
- perr() << "unexpected cell type!\n";
- std::exit(0);
+ perr() << name(descriptor.cell_type_vector[j])
+ << ": unexpected cell type in dimension 3!\n";
+ std::terminate();
+ }
}
}
}
@@ -1672,7 +1670,7 @@ GmshReader::__proceedData()
descriptor.cell_number_vector[jh] = __hexahedra_number[j];
}
- __compute3DCellFaceAndFaceNodeConnectivities(descriptor);
+ __computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
descriptor.cell_owner_vector.resize(nb_cells);
std::fill(descriptor.cell_owner_vector.begin(),
@@ -1691,6 +1689,7 @@ GmshReader::__proceedData()
std::shared_ptr p_connectivity = std::make_shared<Connectivity3D>(descriptor);
Connectivity3D& connectivity = *p_connectivity;
+ const auto& node_number_vector = descriptor.node_number_vector;
using Face = ConnectivityFace<3>;
const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
@@ -1703,7 +1702,7 @@ GmshReader::__proceedData()
for (size_t r=0; r<node_list.size(); ++r) {
node_vector[r] = node_list[r];
}
- face_to_id_map[Face(node_vector)] = l;
+ face_to_id_map[Face(node_vector, node_number_vector)] = l;
}
return face_to_id_map;
} ();
@@ -1712,7 +1711,8 @@ GmshReader::__proceedData()
for (unsigned int f=0; f<__triangles.size(); ++f) {
const unsigned int face_number
= [&]{
- auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}));
+ auto i = face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]},
+ node_number_vector));
if (i == face_to_id_map.end()) {
std::cerr << "face not found!\n";
std::terminate();
@@ -1727,7 +1727,7 @@ GmshReader::__proceedData()
const unsigned int face_number
= [&]{
auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
- __quadrangles[f][2], __quadrangles[f][3]}));
+ __quadrangles[f][2], __quadrangles[f][3]}, node_number_vector));
if (i == face_to_id_map.end()) {
std::cerr << "face not found!\n";
std::terminate();
@@ -1789,7 +1789,7 @@ GmshReader::__proceedData()
descriptor.cell_number_vector[jq] = __quadrangles_number[j];
}
- this->__compute2DCellFaceAndFaceNodeConnectivities(descriptor);
+ this->__computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
descriptor.cell_owner_vector.resize(nb_cells);
std::fill(descriptor.cell_owner_vector.begin(),
@@ -1810,13 +1810,14 @@ GmshReader::__proceedData()
Connectivity2D& connectivity = *p_connectivity;
using Face = ConnectivityFace<2>;
+ const auto& node_number_vector = descriptor.node_number_vector;
const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map
= [&] {
std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
const auto& face_node_list = connectivity.faceToNodeMatrix();
for (FaceId l=0; l<connectivity.numberOfFaces(); ++l) {
const auto& node_list = face_node_list[l];
- face_to_id_map[Face({node_list[0], node_list[1]})] = l;
+ face_to_id_map[Face({node_list[0], node_list[1]}, node_number_vector)] = l;
}
return face_to_id_map;
} ();
@@ -1825,7 +1826,7 @@ GmshReader::__proceedData()
for (unsigned int e=0; e<__edges.size(); ++e) {
const unsigned int edge_number
= [&]{
- auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}));
+ auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}, node_number_vector));
if (i == face_to_id_map.end()) {
std::cerr << "face " << __edges[e][0] << " not found!\n";
std::terminate();
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index cd7d89f23..1ddb6c4db 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -258,8 +258,9 @@ private:
*
*/
void __readPhysicalNames2_2();
- void __compute2DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
- void __compute3DCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+
+ template <size_t Dimension>
+ void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
template <int Dimension>
void _dispatch();
--
GitLab