From 9721b088b5158c770fb1f43089171718947ed9cf Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Tue, 9 Apr 2019 19:13:19 +0200
Subject: [PATCH] Manage edges in 3D
- building/reading of edges
- dispatching in parallel
This implies that [Sub]ItemValue<edge> are also now available
---
src/mesh/Connectivity.cpp | 60 +++-
src/mesh/Connectivity.hpp | 29 +-
src/mesh/ConnectivityDispatcher.cpp | 50 ++-
src/mesh/ConnectivityDispatcher.hpp | 1 +
src/mesh/GmshReader.cpp | 473 ++++++++++++++++++++--------
src/mesh/GmshReader.hpp | 3 +
6 files changed, 465 insertions(+), 151 deletions(-)
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 5c9d59836..3a35c0255 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -54,7 +54,6 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_cell_global_index = cell_global_index;
}
-
{
WeakCellValue<double> inv_cell_nb_nodes(*this);
parallel_for(this->numberOfCells(), PASTIS_LAMBDA(const CellId& j) {
@@ -94,14 +93,13 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_node_is_owned = node_is_owned;
}
+ m_ref_node_list_vector = descriptor.template refItemListVector<ItemType::node>();
+ m_ref_cell_list_vector = descriptor.template refItemListVector<ItemType::cell>();
+
if constexpr (Dimension>1) {
- auto& face_to_node_matrix
- = m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)];
- face_to_node_matrix = descriptor.face_to_node_vector;
+ m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
- 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_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
{
FaceValuePerCell<bool> cell_face_is_reversed(*this);
@@ -111,14 +109,15 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
cell_face_is_reversed(j,lj) = face_cells_vector[lj];
}
}
-
m_cell_face_is_reversed = cell_face_is_reversed;
}
+
{
WeakFaceValue<int> face_number(*this);
face_number = convert_to_array(descriptor.face_number_vector);
m_face_number = face_number;
}
+
{
WeakFaceValue<int> face_owner(*this);
face_owner = convert_to_array(descriptor.face_owner_vector);
@@ -134,10 +133,49 @@ _buildFrom(const ConnectivityDescriptor& descriptor)
m_face_is_owned = face_is_owned;
}
- m_ref_node_list_vector = descriptor.template refItemListVector<ItemType::node>();
- m_ref_edge_list_vector = descriptor.template refItemListVector<ItemType::edge>();
m_ref_face_list_vector = descriptor.template refItemListVector<ItemType::face>();
- m_ref_cell_list_vector = descriptor.template refItemListVector<ItemType::cell>();
+
+ if constexpr (Dimension>2) {
+ m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
+
+ m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
+
+ m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
+
+ {
+ EdgeValuePerFace<bool> face_edge_is_reversed(*this);
+ for (FaceId l=0; l<descriptor.face_edge_is_reversed_vector.size(); ++l) {
+ const auto& edge_faces_vector = descriptor.face_edge_is_reversed_vector[l];
+ for (unsigned short el=0; el<edge_faces_vector.size(); ++el) {
+ face_edge_is_reversed(l,el) = edge_faces_vector[el];
+ }
+ }
+ m_face_edge_is_reversed = face_edge_is_reversed;
+ }
+
+ {
+ WeakEdgeValue<int> edge_number(*this);
+ edge_number = convert_to_array(descriptor.edge_number_vector);
+ m_edge_number = edge_number;
+ }
+
+ {
+ WeakEdgeValue<int> edge_owner(*this);
+ edge_owner = convert_to_array(descriptor.edge_owner_vector);
+ m_edge_owner = edge_owner;
+ }
+
+ {
+ const int rank = parallel::rank();
+ WeakEdgeValue<bool> edge_is_owned(*this);
+ parallel_for(this->numberOfEdges(), PASTIS_LAMBDA(const EdgeId& e) {
+ edge_is_owned[e] = (m_edge_owner[e] == rank);
+ });
+ m_edge_is_owned = edge_is_owned;
+ }
+
+ m_ref_edge_list_vector = descriptor.template refItemListVector<ItemType::edge>();
+ }
}
}
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index 2a27cb380..098675884 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -48,14 +48,13 @@ class ConnectivityDescriptor
public:
std::vector<std::vector<unsigned int>> cell_to_node_vector;
-
-#warning should be set in 2d only
std::vector<std::vector<unsigned int>> cell_to_face_vector;
-
-#warning should be set in 3d only
- // std::vector<std::vector<unsigned int>> cell_to_edge_vector;
+ std::vector<std::vector<unsigned int>> cell_to_edge_vector;
std::vector<std::vector<unsigned int>> face_to_node_vector;
+ std::vector<std::vector<unsigned int>> face_to_edge_vector;
+
+ std::vector<std::vector<unsigned int>> edge_to_node_vector;
template <typename ItemOfItemT>
auto& itemOfItemVector()
@@ -64,19 +63,27 @@ class ConnectivityDescriptor
return cell_to_node_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,FaceOfCell>) {
return cell_to_face_vector;
+ } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfCell>) {
+ return cell_to_edge_vector;
+ } else if constexpr (std::is_same_v<ItemOfItemT,EdgeOfFace>) {
+ return face_to_edge_vector;
} else if constexpr (std::is_same_v<ItemOfItemT,NodeOfFace>) {
return face_to_node_vector;
+ } else if constexpr (std::is_same_v<ItemOfItemT,NodeOfEdge>) {
+ return edge_to_node_vector;
} else {
static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
}
}
std::vector<Array<bool>> cell_face_is_reversed_vector;
+ std::vector<Array<bool>> face_edge_is_reversed_vector;
std::vector<CellType> cell_type_vector;
std::vector<int> cell_number_vector;
std::vector<int> face_number_vector;
+ std::vector<int> edge_number_vector;
std::vector<int> node_number_vector;
template <ItemType item_type>
@@ -86,6 +93,8 @@ class ConnectivityDescriptor
return cell_number_vector;
} else if constexpr (item_type == ItemType::face) {
return face_number_vector;
+ } else if constexpr (item_type == ItemType::edge) {
+ return edge_number_vector;
} else if constexpr (item_type == ItemType::node) {
return node_number_vector;
} else {
@@ -95,6 +104,7 @@ class ConnectivityDescriptor
std::vector<int> cell_owner_vector;
std::vector<int> face_owner_vector;
+ std::vector<int> edge_owner_vector;
std::vector<int> node_owner_vector;
template <ItemType item_type>
@@ -167,7 +177,6 @@ class Connectivity final
WeakCellValue<const int> m_cell_number;
WeakFaceValue<const int> m_face_number;
-#warning check that m_edge_number is filled correctly
WeakEdgeValue<const int> m_edge_number;
WeakNodeValue<const int> m_node_number;
@@ -184,6 +193,7 @@ class Connectivity final
WeakNodeValue<const bool> m_node_is_owned;
WeakFaceValuePerCell<const bool> m_cell_face_is_reversed;
+ WeakEdgeValuePerFace<const bool> m_face_edge_is_reversed;
WeakNodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
WeakEdgeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_edges;
@@ -474,6 +484,13 @@ class Connectivity final
return m_cell_face_is_reversed;
}
+ PASTIS_INLINE
+ const auto& faceEdgeIsReversed() const
+ {
+ static_assert(Dimension>2, "reversed edges makes no sense in dimension 1 or 2");
+ return m_face_edge_is_reversed;
+ }
+
PASTIS_INLINE
const auto& cellLocalNumbersInTheirNodes() const
{
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index 5b9fcdd4e..a21b316bd 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -431,6 +431,8 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
return number_of_item_list_sender;
}();
+ pout() << "- dispatching references for " << itemName(item_type) << '\n';
+
if (number_of_item_list_sender > 0) {
if (number_of_item_list_sender > 1) {
perr() << __FILE__ << ':' << __LINE__ << ": "
@@ -503,6 +505,10 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
return ref_id_list;
} ();
+ for (auto ref_id : ref_id_list) {
+ pout() << " - " << ref_id << '\n';
+ }
+
using block_type = int32_t;
constexpr size_t block_size = sizeof(block_type);
const size_t nb_block = ref_id_list.size()/block_size + (ref_id_list.size()%block_size != 0);
@@ -575,13 +581,41 @@ ConnectivityDispatcher<Dimension>::_buildItemReferenceList()
m_new_descriptor.addRefItemList(RefItemList<item_type>(ref_id_list[i_ref], item_id_array));
}
-
- pout() << __FILE__ << ':' << __LINE__ << ": remains to build lists\n";
}
}
}
}
+template <int Dimension>
+void
+ConnectivityDispatcher<Dimension>::_dispatchEdges()
+{
+ if constexpr (Dimension>2) {
+ this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfCell>();
+ this->_buildSubItemNumbersToRecvByProc<EdgeOfCell>();
+ this->_buildSubItemNumberToIdMap<EdgeOfCell>();
+ this->_buildItemToExchangeLists<ItemType::edge>();
+
+ this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), m_new_descriptor.edge_number_vector);
+
+ this->_buildItemToSubItemDescriptor<EdgeOfCell>();
+
+ this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfEdge>();
+ this->_buildSubItemNumbersToRecvByProc<NodeOfEdge>();
+ this->_buildItemToSubItemDescriptor<NodeOfEdge>();
+
+ this->_buildNumberOfSubItemPerItemToRecvByProc<EdgeOfFace>();
+ this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
+ this->_buildItemToSubItemDescriptor<EdgeOfFace>();
+
+ this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), m_new_descriptor.face_edge_is_reversed_vector);
+
+ this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, m_new_descriptor.edge_owner_vector);
+
+ this->_buildItemReferenceList<ItemType::edge>();
+ }
+}
+
template <int Dimension>
void
ConnectivityDispatcher<Dimension>::_dispatchFaces()
@@ -592,6 +626,10 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
this->_buildSubItemNumberToIdMap<FaceOfCell>();
this->_buildItemToExchangeLists<ItemType::face>();
+ this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfFace>();
+ this->_buildSubItemNumbersToRecvByProc<NodeOfFace>();
+ this->_buildItemToSubItemDescriptor<NodeOfFace>();
+
this->_gatherFrom(m_connectivity.template number<ItemType::face>(), m_new_descriptor.face_number_vector);
this->_buildItemToSubItemDescriptor<FaceOfCell>();
@@ -600,10 +638,6 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, m_new_descriptor.face_owner_vector);
- this->_buildNumberOfSubItemPerItemToRecvByProc<NodeOfFace>();
- this->_buildSubItemNumbersToRecvByProc<NodeOfFace>();
- this->_buildItemToSubItemDescriptor<NodeOfFace>();
-
this->_buildItemReferenceList<ItemType::face>();
}
}
@@ -615,7 +649,7 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
{
this->_buildNewOwner<ItemType::cell>();
this->_buildNewOwner<ItemType::face>();
- // this->_buildNewOwner<ItemType::edge>();
+ this->_buildNewOwner<ItemType::edge>();
this->_buildNewOwner<ItemType::node>();
this->_buildItemToExchangeLists<ItemType::cell>();
@@ -643,6 +677,8 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
this->_dispatchFaces();
+ this->_dispatchEdges();
+
this->_buildItemReferenceList<ItemType::node>();
m_dispatched_connectivity = ConnectivityType::build(m_new_descriptor);
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index d6a97afcb..710f336a9 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -174,6 +174,7 @@ class ConnectivityDispatcher
template <typename ItemOfItemT>
void _buildItemToSubItemDescriptor();
+ void _dispatchEdges();
void _dispatchFaces();
template<typename DataType, ItemType item_type, typename ConnectivityPtr>
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index d1e4e5412..2817d7a3d 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -937,7 +937,7 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
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];
+ const auto& [cell_number, cell_local_face, reversed] = cells_vector[lj];
descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
}
++l;
@@ -978,6 +978,134 @@ GmshReader::__computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& d
}
+template <size_t Dimension>
+void
+GmshReader::__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+{
+ static_assert(Dimension==3, "Invalid dimension to compute face-edge connectivities");
+ using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
+ using Edge = ConnectivityFace<2>;
+
+ const auto& node_number_vector = descriptor.node_number_vector;
+ Array<unsigned short> face_nb_edges(descriptor.face_to_node_vector.size());
+ std::map<Edge, std::vector<FaceEdgeInfo>> edge_faces_map;
+ for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
+ const auto& face_nodes = descriptor.face_to_node_vector[l];
+
+ face_nb_edges[l] = face_nodes.size();
+ for (size_t r=0; r<face_nodes.size()-1; ++r) {
+ Edge e({face_nodes[r], face_nodes[r+1]}, node_number_vector);
+ edge_faces_map[e].emplace_back(std::make_tuple(l, r, e.reversed()));
+ }
+ {
+ Edge e({face_nodes[face_nodes.size()-1], face_nodes[0]}, node_number_vector);
+ edge_faces_map[e].emplace_back(std::make_tuple(l, face_nodes.size()-1, e.reversed()));
+ }
+ }
+
+ std::unordered_map<Edge,EdgeId,typename Edge::Hash> edge_id_map;
+ {
+ descriptor.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
+ for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
+ descriptor.face_to_edge_vector[l].resize(face_nb_edges[l]);
+ }
+ EdgeId e=0;
+ for (const auto& edge_faces_vector : edge_faces_map) {
+ const auto& faces_vector = edge_faces_vector.second;
+ for (unsigned short l=0; l<faces_vector.size(); ++l) {
+ const auto& [face_number, face_local_edge, reversed] = faces_vector[l];
+ descriptor.face_to_edge_vector[face_number][face_local_edge] = e;
+ }
+ edge_id_map[edge_faces_vector.first] = e;
+ ++e;
+ }
+ }
+
+ {
+ descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
+ for (FaceId j=0; j<descriptor.face_edge_is_reversed_vector.size(); ++j) {
+ descriptor.face_edge_is_reversed_vector[j] = Array<bool>(face_nb_edges[j]);
+ }
+ for (const auto& edge_faces_vector : edge_faces_map) {
+ const auto& faces_vector = edge_faces_vector.second;
+ for (unsigned short lj=0; lj<faces_vector.size(); ++lj) {
+ const auto& [face_number, face_local_edge, reversed] = faces_vector[lj];
+ descriptor.face_edge_is_reversed_vector[face_number][face_local_edge] = reversed;
+ }
+ }
+ }
+
+ {
+ descriptor.edge_to_node_vector.resize(edge_faces_map.size());
+ int e=0;
+ for (const auto& edge_info : edge_faces_map) {
+ const Edge& edge = edge_info.first;
+ descriptor.edge_to_node_vector[e] = edge.nodeIdList();
+ ++e;
+ }
+ }
+
+ {
+ // Edge numbers may change if numbers are provided in the file
+ descriptor.edge_number_vector.resize(edge_faces_map.size());
+ for (size_t e=0; e<edge_faces_map.size(); ++e) {
+ descriptor.edge_number_vector[e] = e;
+ }
+ }
+
+ {
+ descriptor.cell_to_node_vector.reserve(descriptor.cell_to_node_vector.size());
+ for (CellId j=0; j<descriptor.cell_to_node_vector.size(); ++j) {
+ const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Tetrahedron: {
+ constexpr int local_edge[6][2] = {{0,1}, {0,2}, {0,3}, {1,2}, {2,3}, {3,1}};
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(6);
+ for (int i_edge=0; i_edge<6; ++i_edge) {
+ const auto e = local_edge[i_edge];
+ Edge edge{{cell_nodes[e[0]],cell_nodes[e[1]]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ std::cerr << "could not find this edge!\n";
+ std::terminate();
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ break;
+ }
+ case CellType::Hexahedron: {
+ constexpr int local_edge[12][2] = {{0,1}, {1,2}, {2,3}, {3,0},
+ {4,5}, {5,6}, {6,7}, {7,4},
+ {0,4}, {1,5}, {2,6}, {3,7}};
+ std::vector<unsigned int> cell_edge_vector;
+ cell_edge_vector.reserve(12);
+ for (int i_edge=0; i_edge<12; ++i_edge) {
+ const auto e = local_edge[i_edge];
+ Edge edge{{cell_nodes[e[0]],cell_nodes[e[1]]}, node_number_vector};
+ auto i = edge_id_map.find(edge);
+ if (i == edge_id_map.end()) {
+ std::cerr << "could not find this edge!\n";
+ std::terminate();
+ }
+ cell_edge_vector.push_back(i->second);
+ }
+ descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ break;
+ }
+ default: {
+ perr() << name(descriptor.cell_type_vector[j])
+ << ": unexpected cell type in dimension 3!\n";
+ std::terminate();
+ }
+ }
+ }
+ }
+}
+
+
void
GmshReader::__proceedData()
{
@@ -1308,126 +1436,217 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
}
- __computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
-
- descriptor.cell_owner_vector.resize(nb_cells);
- std::fill(descriptor.cell_owner_vector.begin(),
- descriptor.cell_owner_vector.end(),
- parallel::rank());
-
- descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
- std::fill(descriptor.face_owner_vector.begin(),
- descriptor.face_owner_vector.end(),
- parallel::rank());
-
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(),
- descriptor.node_owner_vector.end(),
- parallel::rank());
+ this->__computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
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
- = [&] {
- std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
- for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
- const auto& node_vector = descriptor.face_to_node_vector[l];
- face_to_id_map[Face(node_vector, node_number_vector)] = l;
- }
- return face_to_id_map;
- } ();
-
- std::unordered_map<int, FaceId> face_number_id_map
- = [&] {
- std::unordered_map<int, FaceId> face_number_id_map;
- for (size_t l=0; l< descriptor.face_number_vector.size(); ++l) {
- face_number_id_map[descriptor.face_number_vector[l]] = l;
- }
- Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
- return face_number_id_map;
- } ();
-
- std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
- for (unsigned int f=0; f<__triangles.size(); ++f) {
- const unsigned int face_id
- = [&]{
- 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();
+ {
+ using Face = ConnectivityFace<3>;
+ 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;
+ for (FaceId l=0; l<descriptor.face_to_node_vector.size(); ++l) {
+ const auto& node_vector = descriptor.face_to_node_vector[l];
+ face_to_id_map[Face(node_vector, node_number_vector)] = l;
}
- return i->second;
- }();
-
- const unsigned int& ref = __triangles_ref[f];
- ref_faces_map[ref].push_back(face_id);
-
- if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
- if (auto i_face = face_number_id_map.find(__quadrangles_number[f]);
- i_face != face_number_id_map.end()) {
- const int other_face_id = i_face->second;
- std::swap(descriptor.face_number_vector[face_id],
- descriptor.face_number_vector[other_face_id]);
+ return face_to_id_map;
+ } ();
+
+ std::unordered_map<int, FaceId> face_number_id_map
+ = [&] {
+ std::unordered_map<int, FaceId> face_number_id_map;
+ for (size_t l=0; l< descriptor.face_number_vector.size(); ++l) {
+ face_number_id_map[descriptor.face_number_vector[l]] = l;
+ }
+ Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
+ return face_number_id_map;
+ } ();
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
+ for (unsigned int f=0; f<__triangles.size(); ++f) {
+ const unsigned int face_id
+ = [&]{
+ 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();
+ }
+ return i->second;
+ }();
+
+ const unsigned int& ref = __triangles_ref[f];
+ ref_faces_map[ref].push_back(face_id);
+
+ if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
+ if (auto i_face = face_number_id_map.find(__quadrangles_number[f]);
+ i_face != face_number_id_map.end()) {
+ const int other_face_id = i_face->second;
+ std::swap(descriptor.face_number_vector[face_id],
+ descriptor.face_number_vector[other_face_id]);
+
+ face_number_id_map.erase(descriptor.face_number_vector[face_id]);
+ face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
+
+ face_number_id_map[descriptor.face_number_vector[face_id]]=face_id;
+ face_number_id_map[descriptor.face_number_vector[other_face_id]]=other_face_id;
+ } else {
+ face_number_id_map.erase(descriptor.face_number_vector[face_id]);
+ descriptor.face_number_vector[face_id] = __quadrangles_number[f];
+ face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ }
+ }
+ }
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
+ for (unsigned int f=0; f<__quadrangles.size(); ++f) {
+ const unsigned int face_id
+ = [&]{
+ auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
+ __quadrangles[f][2], __quadrangles[f][3]}, node_number_vector));
+ if (i == face_to_id_map.end()) {
+ std::cerr << "face not found!\n";
+ std::terminate();
+ }
+ return i->second;
+ }();
+
+ const unsigned int& ref = __quadrangles_ref[f];
+ ref_faces_map[ref].push_back(face_id);
+
+ if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
+ if (auto i_face = face_number_id_map.find(__quadrangles_number[f]);
+ i_face != face_number_id_map.end()) {
+ const int other_face_id = i_face->second;
+ std::swap(descriptor.face_number_vector[face_id],
+ descriptor.face_number_vector[other_face_id]);
+
+ face_number_id_map.erase(descriptor.face_number_vector[face_id]);
+ face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
+
+ face_number_id_map[descriptor.face_number_vector[face_id]]=face_id;
+ face_number_id_map[descriptor.face_number_vector[other_face_id]]=other_face_id;
+ } else {
+ face_number_id_map.erase(descriptor.face_number_vector[face_id]);
+ descriptor.face_number_vector[face_id] = __quadrangles_number[f];
+ face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ }
+ }
+ }
- face_number_id_map[descriptor.face_number_vector[face_id]]=face_id;
- face_number_id_map[descriptor.face_number_vector[other_face_id]]=other_face_id;
- } else {
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- descriptor.face_number_vector[face_id] = __quadrangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ for (const auto& ref_face_list : ref_faces_map) {
+ Array<FaceId> face_list(ref_face_list.second.size());
+ for (size_t j=0; j<ref_face_list.second.size(); ++j) {
+ face_list[j]=ref_face_list.second[j];
}
+ const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
+ descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
}
}
-
- for (unsigned int f=0; f<__quadrangles.size(); ++f) {
- const unsigned int face_id
- = [&]{
- auto i = face_to_id_map.find(Face({__quadrangles[f][0], __quadrangles[f][1],
- __quadrangles[f][2], __quadrangles[f][3]}, node_number_vector));
- if (i == face_to_id_map.end()) {
- std::cerr << "face not found!\n";
- std::terminate();
+ this->__computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(descriptor);
+
+ {
+ using Edge = ConnectivityFace<2>;
+ const auto& node_number_vector = descriptor.node_number_vector;
+ const std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map
+ = [&] {
+ std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_to_id_map;
+ for (EdgeId l=0; l<descriptor.edge_to_node_vector.size(); ++l) {
+ const auto& node_vector = descriptor.edge_to_node_vector[l];
+ edge_to_id_map[Edge(node_vector, node_number_vector)] = l;
}
- return i->second;
- }();
-
- const unsigned int& ref = __quadrangles_ref[f];
- ref_faces_map[ref].push_back(face_id);
-
- if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
- if (auto i_face = face_number_id_map.find(__quadrangles_number[f]);
- i_face != face_number_id_map.end()) {
- const int other_face_id = i_face->second;
- std::swap(descriptor.face_number_vector[face_id],
- descriptor.face_number_vector[other_face_id]);
-
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
+ return edge_to_id_map;
+ } ();
+
+ std::unordered_map<int, EdgeId> edge_number_id_map
+ = [&] {
+ std::unordered_map<int, EdgeId> edge_number_id_map;
+ for (size_t l=0; l< descriptor.edge_number_vector.size(); ++l) {
+ edge_number_id_map[descriptor.edge_number_vector[l]] = l;
+ }
+ Assert(edge_number_id_map.size() == descriptor.edge_number_vector.size());
+ return edge_number_id_map;
+ } ();
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_edges_map;
+ for (unsigned int e=0; e<__edges.size(); ++e) {
+ const unsigned int edge_id
+ = [&]{
+ auto i = edge_to_id_map.find(Edge({__edges[e][0], __edges[e][1]}, node_number_vector));
+ if (i == edge_to_id_map.end()) {
+ std::cerr << "edge " << __edges[e][0] << " not found!\n";
+ std::terminate();
+ }
+ return i->second;
+ }();
+ const unsigned int& ref = __edges_ref[e];
+ ref_edges_map[ref].push_back(edge_id);
+
+ if (descriptor.edge_number_vector[edge_id] != __edges_number[e]) {
+ if (auto i_edge = edge_number_id_map.find(__edges_number[e]);
+ i_edge != edge_number_id_map.end()) {
+ const int other_edge_id = i_edge->second;
+ std::swap(descriptor.edge_number_vector[edge_id],
+ descriptor.edge_number_vector[other_edge_id]);
+
+ edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
+ edge_number_id_map.erase(descriptor.edge_number_vector[other_edge_id]);
+
+ edge_number_id_map[descriptor.edge_number_vector[edge_id]]=edge_id;
+ edge_number_id_map[descriptor.edge_number_vector[other_edge_id]]=other_edge_id;
+ } else {
+ edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
+ descriptor.edge_number_vector[edge_id] = __edges_number[e];
+ edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_id;
+ }
+ }
+ }
- face_number_id_map[descriptor.face_number_vector[face_id]]=face_id;
- face_number_id_map[descriptor.face_number_vector[other_face_id]]=other_face_id;
- } else {
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- descriptor.face_number_vector[face_id] = __quadrangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ for (const auto& ref_edge_list : ref_edges_map) {
+ Array<EdgeId> edge_list(ref_edge_list.second.size());
+ for (size_t j=0; j<ref_edge_list.second.size(); ++j) {
+ edge_list[j]=ref_edge_list.second[j];
}
+ const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_edge_list.first);
+ descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
}
}
- for (const auto& ref_face_list : ref_faces_map) {
- Array<FaceId> face_list(ref_face_list.second.size());
- for (size_t j=0; j<ref_face_list.second.size(); ++j) {
- face_list[j]=ref_face_list.second[j];
+ std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+ for (unsigned int r=0; r<__points.size(); ++r) {
+ const unsigned int point_number = __points[r];
+ const unsigned int& ref = __points_ref[r];
+ ref_points_map[ref].push_back(point_number);
+ }
+
+ for (const auto& ref_point_list : ref_points_map) {
+ Array<NodeId> point_list(ref_point_list.second.size());
+ for (size_t j=0; j<ref_point_list.second.size(); ++j) {
+ point_list[j]=ref_point_list.second[j];
}
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
- descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
+ const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_point_list.first);
+ descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list));
}
+ descriptor.cell_owner_vector.resize(nb_cells);
+ std::fill(descriptor.cell_owner_vector.begin(),
+ descriptor.cell_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
+ std::fill(descriptor.face_owner_vector.begin(),
+ descriptor.face_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.edge_owner_vector.resize(descriptor.edge_number_vector.size());
+ std::fill(descriptor.edge_owner_vector.begin(),
+ descriptor.edge_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(),
+ descriptor.node_owner_vector.end(),
+ parallel::rank());
+
std::shared_ptr p_connectivity = Connectivity3D::build(descriptor);
Connectivity3D& connectivity = *p_connectivity;
@@ -1497,21 +1716,6 @@ GmshReader::__proceedData()
this->__computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
- descriptor.cell_owner_vector.resize(nb_cells);
- std::fill(descriptor.cell_owner_vector.begin(),
- descriptor.cell_owner_vector.end(),
- parallel::rank());
-
- descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
- std::fill(descriptor.face_owner_vector.begin(),
- descriptor.face_owner_vector.end(),
- parallel::rank());
-
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(),
- descriptor.node_owner_vector.end(),
- parallel::rank());
-
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
@@ -1593,6 +1797,20 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list));
}
+ descriptor.cell_owner_vector.resize(nb_cells);
+ std::fill(descriptor.cell_owner_vector.begin(),
+ descriptor.cell_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.face_owner_vector.resize(descriptor.face_number_vector.size());
+ std::fill(descriptor.face_owner_vector.begin(),
+ descriptor.face_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(),
+ descriptor.node_owner_vector.end(),
+ parallel::rank());
std::shared_ptr p_connectivity = Connectivity2D::build(descriptor);
Connectivity2D& connectivity = *p_connectivity;
@@ -1626,15 +1844,6 @@ GmshReader::__proceedData()
descriptor.cell_number_vector[j] = __edges_number[j];
}
- descriptor.cell_owner_vector.resize(nb_cells);
- std::fill(descriptor.cell_owner_vector.begin(),
- descriptor.cell_owner_vector.end(),
- parallel::rank());
-
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(),
- descriptor.node_owner_vector.end(),
- parallel::rank());
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r=0; r<__points.size(); ++r) {
@@ -1652,6 +1861,16 @@ GmshReader::__proceedData()
descriptor.addRefItemList(RefNodeList(physical_ref_id.refId(), point_list));
}
+ descriptor.cell_owner_vector.resize(nb_cells);
+ std::fill(descriptor.cell_owner_vector.begin(),
+ descriptor.cell_owner_vector.end(),
+ parallel::rank());
+
+ descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(),
+ descriptor.node_owner_vector.end(),
+ parallel::rank());
+
std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
Connectivity1D& connectivity = *p_connectivity;
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 1ddb6c4db..cd482e637 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -262,6 +262,9 @@ private:
template <size_t Dimension>
void __computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
+ template <size_t Dimension>
+ void __computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
+
template <int Dimension>
void _dispatch();
--
GitLab