diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 7c7f7c72669a95bd096c697c8aa1929f134dbcda..4ec8145a85513bb68d452f1ffcb28679932ea5dc 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 fda25db04737cab3e10d6047870e6ebc377d11c2..68e4dd3c467594eae95897b27cfac74ea99cc5da 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 cd7d89f236272cdf6c584a270587677c1339f379..1ddb6c4db32dbb3deac2fe6f87364532f4da84aa 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();