diff --git a/doc/userdoc.org b/doc/userdoc.org
index c338582d59ef05c6b2241c09b615a515e56b2fe0..ba7cd798a6a7991d6451c405a5b9c81dfa2c0362 100644
--- a/doc/userdoc.org
+++ b/doc/userdoc.org
@@ -1223,8 +1223,7 @@ are sorted by type of left hand side variable.
Observe that for these small matrix types ($\mathbb{R}^{d\times d}$) the
construction ~A *= B;~ where ~B~ is a matrix of the same type as ~A~ is not
allowed. The main reason for that is that for $d>1$ this operation has
-no interest since it requires a temporary. One will see bellow that it
-is possible to write ~A = A*B;~ if needed.
+no interest since it requires a temporary. One will see below that it
#+END_note
- ~string~: the ~*=~ operator is not defined for left hand side string variables.
diff --git a/src/language/modules/MathFunctionRegisterForVh.cpp b/src/language/modules/MathFunctionRegisterForVh.cpp
index 122f9d20d1eb21c34bbb57764789880ec632583d..38dc08b6e4bec82be8cab2666080563bd7cb9736 100644
--- a/src/language/modules/MathFunctionRegisterForVh.cpp
+++ b/src/language/modules/MathFunctionRegisterForVh.cpp
@@ -227,6 +227,12 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
));
+ scheme_module._addBuiltinFunction("trace", std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> A)
+ -> std::shared_ptr<const DiscreteFunctionVariant> { return trace(A); }
+
+ ));
scheme_module._addBuiltinFunction("min",
std::function(
@@ -234,6 +240,14 @@ MathFunctionRegisterForVh::MathFunctionRegisterForVh(SchemeModule& scheme_module
));
+ scheme_module._addBuiltinFunction("transpose",
+ std::function(
+
+ [](std::shared_ptr<const DiscreteFunctionVariant> A)
+ -> std::shared_ptr<const DiscreteFunctionVariant> { return transpose(A); }
+
+ ));
+
scheme_module._addBuiltinFunction("min", std::function(
[](std::shared_ptr<const DiscreteFunctionVariant> a,
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index 9c5cf3aec1e11d16f652357db53f68db54f6ed79..af1fca45e7dec5898a12731430358691728c7a2a 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -12,6 +12,7 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
#include <mesh/DualMeshManager.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/IBoundaryDescriptor.hpp>
@@ -143,6 +144,37 @@ MeshModule::MeshModule()
));
+ this->_addBuiltinFunction("check_connectivity_ordering",
+ std::function(
+
+ [](const std::shared_ptr<const IMesh>& i_mesh) -> bool {
+ switch (i_mesh->dimension()) {
+ case 1: {
+ using MeshType = Mesh<Connectivity<1>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return checkConnectivityOrdering(p_mesh->connectivity());
+ }
+ case 2: {
+ using MeshType = Mesh<Connectivity<2>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return checkConnectivityOrdering(p_mesh->connectivity());
+ }
+ case 3: {
+ using MeshType = Mesh<Connectivity<3>>;
+
+ std::shared_ptr p_mesh = std::dynamic_pointer_cast<const MeshType>(i_mesh);
+ return checkConnectivityOrdering(p_mesh->connectivity());
+ }
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ }
+ }
+
+ ));
+
this->_addBuiltinFunction("cartesianMesh",
std::function(
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 66ce42ae6f9e490a0b3e6e406934478b0f026f4a..62dbd50536ad2e725e67a25d9cddf481a62b40c6 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -7,6 +7,7 @@ add_library(
ConnectivityBuilderBase.cpp
ConnectivityComputer.cpp
ConnectivityDispatcher.cpp
+ ConnectivityUtils.cpp
DiamondDualConnectivityBuilder.cpp
DiamondDualMeshBuilder.cpp
Dual1DConnectivityBuilder.cpp
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index 22e722a11f0a53aabbde7a243386f854fcd41a59..dc03c1ee65224a07466e5f08352b53004ad3c1ea 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -13,44 +13,36 @@ template <size_t Dimension>
void
Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
{
- Assert(descriptor.cell_to_node_vector.size() == descriptor.cell_type_vector.size());
- Assert(descriptor.cell_number_vector.size() == descriptor.cell_type_vector.size());
+ Assert(descriptor.cellToNodeMatrix().numberOfRows() == descriptor.cellTypeVector().size());
+ Assert(descriptor.cellNumberVector().size() == descriptor.cellTypeVector().size());
if constexpr (Dimension > 1) {
- Assert(descriptor.cell_to_face_vector.size() == descriptor.cell_type_vector.size());
- Assert(descriptor.face_to_node_vector.size() == descriptor.face_number_vector.size());
- Assert(descriptor.face_owner_vector.size() == descriptor.face_number_vector.size());
+ Assert(descriptor.cellToFaceMatrix().numberOfRows() == descriptor.cellTypeVector().size());
+ Assert(descriptor.faceToNodeMatrix().numberOfRows() == descriptor.faceNumberVector().size());
+ Assert(descriptor.faceOwnerVector().size() == descriptor.faceNumberVector().size());
}
- m_number_of_cells = descriptor.cell_number_vector.size();
- m_number_of_nodes = descriptor.node_number_vector.size();
+ m_number_of_cells = descriptor.cellNumberVector().size();
+ m_number_of_nodes = descriptor.nodeNumberVector().size();
if constexpr (Dimension == 1) {
m_number_of_edges = m_number_of_nodes;
m_number_of_faces = m_number_of_nodes;
} else {
- m_number_of_faces = descriptor.face_number_vector.size();
+ m_number_of_faces = descriptor.faceNumberVector().size();
if constexpr (Dimension == 2) {
m_number_of_edges = m_number_of_faces;
} else {
static_assert(Dimension == 3, "unexpected dimension");
- m_number_of_edges = descriptor.edge_number_vector.size();
+ m_number_of_edges = descriptor.edgeNumberVector().size();
}
}
auto& cell_to_node_matrix = m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::node)];
- cell_to_node_matrix = descriptor.cell_to_node_vector;
+ cell_to_node_matrix = descriptor.cellToNodeMatrix();
- {
- WeakCellValue<CellType> cell_type(*this);
- parallel_for(
- this->numberOfCells(), PUGS_LAMBDA(CellId j) { cell_type[j] = descriptor.cell_type_vector[j]; });
- m_cell_type = cell_type;
- }
-
- m_cell_number = WeakCellValue<int>(*this, convert_to_array(descriptor.cell_number_vector));
-
- Array node_number_array = convert_to_array(descriptor.node_number_vector);
- m_node_number = WeakNodeValue<int>(*this, node_number_array);
+ m_cell_type = WeakCellValue<const CellType>(*this, descriptor.cellTypeVector());
+ m_cell_number = WeakCellValue<const int>(*this, descriptor.cellNumberVector());
+ m_node_number = WeakNodeValue<const int>(*this, descriptor.nodeNumberVector());
{
WeakCellValue<int> cell_global_index(*this);
@@ -60,7 +52,7 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
m_cell_global_index = cell_global_index;
}
- m_cell_owner = WeakCellValue<int>(*this, convert_to_array(descriptor.cell_owner_vector));
+ m_cell_owner = WeakCellValue<const int>(*this, descriptor.cellOwnerVector());
{
const int rank = parallel::rank();
@@ -70,8 +62,7 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
m_cell_is_owned = cell_is_owned;
}
- Array node_owner_array = convert_to_array(descriptor.node_owner_vector);
- m_node_owner = WeakNodeValue<int>{*this, node_owner_array};
+ m_node_owner = WeakNodeValue<const int>{*this, descriptor.nodeOwnerVector()};
Array<bool> node_is_owned_array(this->numberOfNodes());
{
@@ -87,13 +78,13 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
if constexpr (Dimension == 1) {
// faces are similar to nodes
- m_face_number = WeakFaceValue<int>(*this, node_number_array);
- m_face_owner = WeakFaceValue<int>(*this, node_owner_array);
+ m_face_number = WeakFaceValue<const int>(*this, descriptor.nodeNumberVector());
+ m_face_owner = WeakFaceValue<const int>(*this, descriptor.nodeOwnerVector());
m_face_is_owned = WeakFaceValue<bool>(*this, node_is_owned_array);
// edges are similar to nodes
- m_edge_number = WeakEdgeValue<int>(*this, node_number_array);
- m_edge_owner = WeakEdgeValue<int>(*this, node_owner_array);
+ m_edge_number = WeakEdgeValue<const int>(*this, descriptor.nodeNumberVector());
+ m_edge_owner = WeakEdgeValue<const int>(*this, descriptor.nodeOwnerVector());
m_edge_is_owned = WeakEdgeValue<bool>(*this, node_is_owned_array);
// edge and face references are set equal to node references
@@ -114,26 +105,15 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
}
} else {
- m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.face_to_node_vector;
+ m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::node)] = descriptor.faceToNodeMatrix();
- m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cell_to_face_vector;
+ m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::face)] = descriptor.cellToFaceMatrix();
- {
- 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;
- }
+ m_cell_face_is_reversed = FaceValuePerCell<const bool>(*this, descriptor.cellFaceIsReversed());
- Array face_number_array = convert_to_array(descriptor.face_number_vector);
- m_face_number = WeakFaceValue<int>(*this, face_number_array);
+ m_face_number = WeakFaceValue<const int>(*this, descriptor.faceNumberVector());
- Array face_owner_array = convert_to_array(descriptor.face_owner_vector);
- m_face_owner = WeakFaceValue<int>(*this, face_owner_array);
+ m_face_owner = WeakFaceValue<const int>(*this, descriptor.faceOwnerVector());
Array<bool> face_is_owned_array(this->numberOfFaces());
{
@@ -148,8 +128,8 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
if constexpr (Dimension == 2) {
// edges are similar to faces
- m_edge_number = WeakEdgeValue<int>(*this, face_number_array);
- m_edge_owner = WeakEdgeValue<int>(*this, face_owner_array);
+ m_edge_number = WeakEdgeValue<const int>(*this, descriptor.faceNumberVector());
+ m_edge_owner = WeakEdgeValue<const int>(*this, descriptor.faceOwnerVector());
m_edge_is_owned = WeakEdgeValue<bool>(*this, face_is_owned_array);
// edge references are set equal to face references
@@ -166,25 +146,16 @@ Connectivity<Dimension>::_buildFrom(const ConnectivityDescriptor& descriptor)
}
} else {
- m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edge_to_node_vector;
+ m_item_to_item_matrix[itemTId(ItemType::edge)][itemTId(ItemType::node)] = descriptor.edgeToNodeMatrix();
- m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.face_to_edge_vector;
+ m_item_to_item_matrix[itemTId(ItemType::face)][itemTId(ItemType::edge)] = descriptor.faceToEdgeMatrix();
- m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cell_to_edge_vector;
+ m_item_to_item_matrix[itemTId(ItemType::cell)][itemTId(ItemType::edge)] = descriptor.cellToEdgeMatrix();
- {
- 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;
- }
+ m_face_edge_is_reversed = EdgeValuePerFace<const bool>(*this, descriptor.faceEdgeIsReversed());
- m_edge_number = WeakEdgeValue<int>(*this, convert_to_array(descriptor.edge_number_vector));
- m_edge_owner = WeakEdgeValue<int>(*this, convert_to_array(descriptor.edge_owner_vector));
+ m_edge_number = WeakEdgeValue<const int>(*this, descriptor.edgeNumberVector());
+ m_edge_owner = WeakEdgeValue<const int>(*this, descriptor.edgeOwnerVector());
{
const int rank = parallel::rank();
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index b0279507ee41e6f6d3433bda2b5dd7beaf406265..ee7b07fe043a44262c23afc053e748bcac448b50 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -122,7 +122,7 @@ class Connectivity final : public IConnectivity
const ConnectivityMatrix& connectivity_matrix = m_item_to_item_matrix[itemTId(item_type_0)][itemTId(item_type_1)];
if (not connectivity_matrix.isBuilt()) {
const_cast<ConnectivityMatrix&>(connectivity_matrix) =
- m_connectivity_computer.computeConnectivityMatrix(*this, item_type_0, item_type_1);
+ m_connectivity_computer.computeInverseConnectivityMatrix(*this, item_type_0, item_type_1);
}
return connectivity_matrix;
}
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
index 55c56b3a4b78f28ecdf8613d04d4d0348b9e5c4a..e7eb95fc6ed064805d12b5a8ece6e2844c37b9e9 100644
--- a/src/mesh/ConnectivityBuilderBase.cpp
+++ b/src/mesh/ConnectivityBuilderBase.cpp
@@ -7,8 +7,6 @@
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
-#include <map>
-#include <unordered_map>
#include <vector>
template <size_t Dimension>
@@ -16,240 +14,583 @@ void
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
{
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>;
- const auto& node_number_vector = descriptor.node_number_vector;
- Array<unsigned short> cell_nb_faces(descriptor.cell_to_node_vector.size());
- std::map<Face, std::vector<CellFaceInfo>> face_cells_map;
- for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
- const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+ const auto& node_number_vector = descriptor.nodeNumberVector();
+ const auto& cell_to_node_matrix = descriptor.cellToNodeMatrix();
+ const auto& cell_type_vector = descriptor.cellTypeVector();
+
+ size_t total_number_of_faces = 0;
+
+ for (CellId j = 0; j < cell_to_node_matrix.numberOfRows(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
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()));
-
- // 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()));
-
- // 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()));
+ total_number_of_faces += cell_nodes.size();
+ } else if constexpr (Dimension == 3) {
+ switch (cell_type_vector[j]) {
+ case CellType::Hexahedron: {
+ total_number_of_faces += 6;
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()));
-
- // 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()));
-
- // 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()));
-
- // 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()));
+ case CellType::Tetrahedron: {
+ total_number_of_faces += 4;
break;
}
- case CellType::Polygon: {
- cell_nb_faces[j] = cell_nodes.size();
- for (size_t i = 0; i < cell_nodes.size(); ++i) {
- Face f({cell_nodes[i], cell_nodes[(i + 1) % cell_nodes.size()]}, node_number_vector);
- face_cells_map[f].emplace_back(std::make_tuple(j, i, f.reversed()));
- }
+ case CellType::Prism: {
+ total_number_of_faces += 5;
+ break;
+ }
+ case CellType::Pyramid: {
+ total_number_of_faces += cell_nodes.size();
+ break;
+ }
+ case CellType::Diamond: {
+ total_number_of_faces += 2 * (cell_nodes.size() - 2);
break;
}
default: {
std::ostringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+ error_msg << name(cell_type_vector[j]) << ": unexpected cell type in dimension 3";
throw UnexpectedError(error_msg.str());
}
}
- } else if constexpr (Dimension == 3) {
- switch (descriptor.cell_type_vector[j]) {
- case CellType::Hexahedron: {
- // face 0
- Face f0({cell_nodes[3], cell_nodes[2], cell_nodes[1], 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]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+ const size_t total_number_of_face_by_node = [&] {
+ if constexpr (Dimension == 2) {
+ return 2 * total_number_of_faces;
+ } else {
+ Assert(Dimension == 3);
+ size_t count_number_of_face_by_node = 0;
+ for (CellId j = 0; j < cell_to_node_matrix.numberOfRows(); ++j) {
+ switch (cell_type_vector[j]) {
+ case CellType::Hexahedron: {
+ count_number_of_face_by_node += 6 * 4;
+ break;
+ }
+ case CellType::Tetrahedron: {
+ count_number_of_face_by_node += 4 * 3;
+ break;
+ }
+ case CellType::Prism: {
+ count_number_of_face_by_node += 3 * 4 + 2 * 3;
+ break;
+ }
+ case CellType::Pyramid: {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ count_number_of_face_by_node += 1 * cell_nodes.size() + cell_nodes.size() * 3;
+ break;
+ }
+ case CellType::Diamond: {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+ count_number_of_face_by_node += cell_nodes.size() * 3 * 2;
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << name(cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ }
+ return count_number_of_face_by_node;
+ }
+ }();
- // face 2
- Face f2({cell_nodes[0], cell_nodes[4], cell_nodes[7], cell_nodes[3]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+ Array<unsigned int> dup_faces_to_node_list(total_number_of_face_by_node);
+ Array<unsigned int> dup_face_to_node_row(total_number_of_faces + 1);
+ size_t i_face = 0;
+ dup_face_to_node_row[0] = 0;
- // face 3
- Face f3({cell_nodes[1], cell_nodes[2], cell_nodes[6], cell_nodes[5]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ Array<unsigned short> cell_nb_faces(cell_to_node_matrix.numberOfRows());
+ {
+ size_t i_face_node = 0;
+ for (CellId j = 0; j < cell_to_node_matrix.numberOfRows(); ++j) {
+ const auto& cell_nodes = cell_to_node_matrix[j];
+
+ if constexpr (Dimension == 2) {
+ switch (cell_type_vector[j]) {
+ case CellType::Triangle: {
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ cell_nb_faces[j] = 3;
+ break;
+ }
+ case CellType::Quadrangle: {
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ cell_nb_faces[j] = 4;
+ break;
+ }
+ case CellType::Polygon: {
+ for (size_t i = 0; i < cell_nodes.size(); ++i) {
+ dup_faces_to_node_list[i_face_node] = cell_nodes[i];
+ i_face_node++;
+ dup_faces_to_node_list[i_face_node] = cell_nodes[(i + 1) % cell_nodes.size()];
+ i_face_node++;
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 2;
+ i_face++;
+ }
- // face 4
- Face f4({cell_nodes[0], cell_nodes[1], cell_nodes[5], cell_nodes[4]}, node_number_vector);
- face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+ cell_nb_faces[j] = cell_nodes.size();
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << name(cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ } else if constexpr (Dimension == 3) {
+ switch (cell_type_vector[j]) {
+ case CellType::Hexahedron: {
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[7];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[6];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ cell_nb_faces[j] = 6;
+ break;
+ }
+ case CellType::Tetrahedron: {
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
- // face 5
- Face f5({cell_nodes[3], cell_nodes[7], cell_nodes[6], cell_nodes[2]}, node_number_vector);
- face_cells_map[f5].emplace_back(std::make_tuple(j, 5, f5.reversed()));
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
- cell_nb_faces[j] = 6;
- break;
- }
- case CellType::Tetrahedron: {
- cell_nb_faces[j] = 4;
- // face 0
- Face f0({cell_nodes[1], cell_nodes[2], 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]}, 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]}, 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]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
- break;
- }
- case CellType::Prism: {
- // face 0
- Face f0({cell_nodes[2], cell_nodes[1], cell_nodes[0]}, node_number_vector);
- face_cells_map[f0].emplace_back(std::make_tuple(j, 0, f0.reversed()));
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
- // face 1
- Face f1({cell_nodes[3], cell_nodes[4], cell_nodes[5]}, node_number_vector);
- face_cells_map[f1].emplace_back(std::make_tuple(j, 1, f1.reversed()));
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
- // face 2
- Face f2({cell_nodes[1], cell_nodes[2], cell_nodes[5], cell_nodes[4]}, node_number_vector);
- face_cells_map[f2].emplace_back(std::make_tuple(j, 2, f2.reversed()));
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
- // face 3
- Face f3({cell_nodes[0], cell_nodes[1], cell_nodes[4], cell_nodes[3]}, node_number_vector);
- face_cells_map[f3].emplace_back(std::make_tuple(j, 3, f3.reversed()));
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
- // face 4
- Face f4({cell_nodes[2], cell_nodes[0], cell_nodes[3], cell_nodes[5]}, node_number_vector);
- face_cells_map[f4].emplace_back(std::make_tuple(j, 4, f4.reversed()));
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
- cell_nb_faces[j] = 5;
- break;
- }
- case CellType::Pyramid: {
- cell_nb_faces[j] = cell_nodes.size();
- std::vector<unsigned int> base_nodes(cell_nodes.size() - 1);
- for (size_t i = 0; i < base_nodes.size(); ++i) {
- base_nodes[base_nodes.size() - 1 - i] = cell_nodes[i];
- }
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
- // base face
- {
- Face base_face(base_nodes, node_number_vector);
- face_cells_map[base_face].emplace_back(std::make_tuple(j, 0, base_face.reversed()));
+ cell_nb_faces[j] = 4;
+ break;
}
- // side faces
- const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
- for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
- Face side_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], pyramid_vertex},
- node_number_vector);
- face_cells_map[side_face].emplace_back(std::make_tuple(j, i_node + 1, side_face.reversed()));
+ case CellType::Prism: {
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[1];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[4];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[2];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[0];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[3];
+ dup_faces_to_node_list[i_face_node++] = cell_nodes[5];
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 4;
+ i_face++;
+
+ cell_nb_faces[j] = 5;
+ break;
}
- break;
- }
- case CellType::Diamond: {
- cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
- std::vector<unsigned int> base_nodes;
- std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+ case CellType::Pyramid: {
+ cell_nb_faces[j] = cell_nodes.size();
+ std::vector<unsigned int> base_nodes(cell_nodes.size() - 1);
+ for (size_t i = 0; i < base_nodes.size(); ++i) {
+ base_nodes[base_nodes.size() - 1 - i] = cell_nodes[i];
+ }
- { // top faces
- const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i = 0; i < base_nodes.size(); ++i) {
+ dup_faces_to_node_list[i_face_node++] = base_nodes[i];
+ }
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + base_nodes.size();
+ i_face++;
+
+ // side faces
+ const auto pyramid_vertex = cell_nodes[cell_nodes.size() - 1];
for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
- Face top_face({base_nodes[i_node], base_nodes[(i_node + 1) % base_nodes.size()], top_vertex},
- node_number_vector);
- face_cells_map[top_face].emplace_back(std::make_tuple(j, i_node, top_face.reversed()));
+ dup_faces_to_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
+ dup_faces_to_node_list[i_face_node++] = base_nodes[i_node];
+ dup_faces_to_node_list[i_face_node++] = pyramid_vertex;
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
}
+ break;
}
+ case CellType::Diamond: {
+ auto base_nodes = [&](size_t i) { return cell_nodes[i + 1]; };
+
+ { // top faces
+ const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_node = 0; i_node < cell_nodes.size() - 2; ++i_node) {
+ dup_faces_to_node_list[i_face_node++] = base_nodes(i_node);
+ dup_faces_to_node_list[i_face_node++] = base_nodes((i_node + 1) % (cell_nodes.size() - 2));
+ dup_faces_to_node_list[i_face_node++] = top_vertex;
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
+ }
+ }
- { // bottom faces
- const auto bottom_vertex = cell_nodes[0];
- for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
- Face bottom_face({base_nodes[(i_node + 1) % base_nodes.size()], base_nodes[i_node], bottom_vertex},
- node_number_vector);
- face_cells_map[bottom_face].emplace_back(
- std::make_tuple(j, i_node + base_nodes.size(), bottom_face.reversed()));
+ { // bottom faces
+ const auto bottom_vertex = cell_nodes[0];
+ for (size_t i_node = 0; i_node < cell_nodes.size() - 2; ++i_node) {
+ dup_faces_to_node_list[i_face_node++] = base_nodes((i_node + 1) % (cell_nodes.size() - 2));
+ dup_faces_to_node_list[i_face_node++] = base_nodes(i_node);
+ dup_faces_to_node_list[i_face_node++] = bottom_vertex;
+
+ dup_face_to_node_row[i_face + 1] = dup_face_to_node_row[i_face] + 3;
+ i_face++;
+ }
}
+ cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
+ break;
+ }
+ default: {
+ std::ostringstream error_msg;
+ error_msg << name(cell_type_vector[j]) << ": unexpected cell type in dimension 3";
+ throw UnexpectedError(error_msg.str());
+ }
}
- break;
}
- default: {
- std::ostringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
- throw UnexpectedError(error_msg.str());
+ }
+ }
+
+ Array<bool> cell_face_is_reversed(total_number_of_faces);
+
+ if constexpr (Dimension == 2) {
+ for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
+ if (node_number_vector[dup_faces_to_node_list[2 * i_face]] >
+ node_number_vector[dup_faces_to_node_list[2 * i_face + 1]]) {
+ std::swap(dup_faces_to_node_list[2 * i_face], dup_faces_to_node_list[2 * i_face + 1]);
+ cell_face_is_reversed[i_face] = true;
+ } else {
+ cell_face_is_reversed[i_face] = false;
}
+ }
+ } else if constexpr (Dimension == 3) {
+ std::vector<int> buffer;
+
+ for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
+ const size_t face_node_number = dup_face_to_node_row[i_face + 1] - dup_face_to_node_row[i_face];
+ size_t i_face_node_smallest_number = 0;
+ for (size_t i_face_node = 1; i_face_node < face_node_number; ++i_face_node) {
+ if (node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_face_node]] <
+ node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_face_node_smallest_number]]) {
+ i_face_node_smallest_number = i_face_node;
+ }
}
+
+ if (i_face_node_smallest_number != 0) {
+ buffer.resize(face_node_number);
+ for (size_t i_node = i_face_node_smallest_number; i_node < face_node_number; ++i_node) {
+ buffer[i_node - i_face_node_smallest_number] = dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node];
+ }
+ for (size_t i_node = 0; i_node < i_face_node_smallest_number; ++i_node) {
+ buffer[i_node + face_node_number - i_face_node_smallest_number] =
+ dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node];
+ }
+
+ for (size_t i_node = 0; i_node < face_node_number; ++i_node) {
+ dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node] = buffer[i_node];
+ }
+ }
+
+ if (node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face] + 1]] >
+ node_number_vector[dup_faces_to_node_list[dup_face_to_node_row[i_face + 1] - 1]]) {
+ for (size_t i_node = 1; i_node <= (face_node_number + 1) / 2 - 1; ++i_node) {
+ std::swap(dup_faces_to_node_list[dup_face_to_node_row[i_face] + i_node],
+ dup_faces_to_node_list[dup_face_to_node_row[i_face + 1] - i_node]);
+ }
+
+ cell_face_is_reversed[i_face] = true;
+ } else {
+ cell_face_is_reversed[i_face] = false;
+ }
+ }
+ }
+
+ Array<unsigned int> node_to_duplicate_face_row(node_number_vector.size() + 1);
+ node_to_duplicate_face_row.fill(0);
+ for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
+ for (size_t i_node_face = dup_face_to_node_row[i_face]; i_node_face < dup_face_to_node_row[i_face + 1];
+ ++i_node_face) {
+ node_to_duplicate_face_row[dup_faces_to_node_list[i_node_face] + 1]++;
}
}
+ for (size_t i_node = 1; i_node < node_to_duplicate_face_row.size(); ++i_node) {
+ node_to_duplicate_face_row[i_node] += node_to_duplicate_face_row[i_node - 1];
+ }
+
+ Array<unsigned int> node_duplicate_face_list(node_to_duplicate_face_row[node_to_duplicate_face_row.size() - 1]);
+
{
- descriptor.cell_to_face_vector.resize(descriptor.cell_to_node_vector.size());
- for (CellId j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
- descriptor.cell_to_face_vector[j].resize(cell_nb_faces[j]);
+ Array<unsigned int> node_duplicate_face_row_idx(node_number_vector.size());
+ node_duplicate_face_row_idx.fill(0);
+
+ for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
+ for (size_t i_node_face = dup_face_to_node_row[i_face]; i_node_face < dup_face_to_node_row[i_face + 1];
+ ++i_node_face) {
+ const size_t node_id = dup_faces_to_node_list[i_node_face];
+ node_duplicate_face_list[node_to_duplicate_face_row[node_id] + node_duplicate_face_row_idx[node_id]] = i_face;
+ node_duplicate_face_row_idx[node_id]++;
+ }
}
- FaceId 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];
- descriptor.cell_to_face_vector[cell_number][cell_local_face] = l;
- }
- ++l;
+ }
+
+ Array<unsigned int> dup_face_to_face(total_number_of_faces);
+ parallel_for(
+ total_number_of_faces, PUGS_LAMBDA(size_t i_face) { dup_face_to_face[i_face] = i_face; });
+
+ auto is_same_face = [=](const size_t i_face, const size_t j_face) {
+ if ((dup_face_to_node_row[i_face + 1] - dup_face_to_node_row[i_face]) !=
+ (dup_face_to_node_row[j_face + 1] - dup_face_to_node_row[j_face])) {
+ return false;
+ } else {
+ auto i_face_node = dup_face_to_node_row[i_face];
+ auto j_face_node = dup_face_to_node_row[j_face];
+ while (i_face_node < dup_face_to_node_row[i_face + 1]) {
+ if (dup_faces_to_node_list[i_face_node] != dup_faces_to_node_list[j_face_node]) {
+ return false;
+ }
+ i_face_node++;
+ j_face_node++;
+ }
+ return true;
+ }
+ };
+
+ size_t nb_dup_faces = 0;
+ for (size_t i_node = 0; i_node < node_number_vector.size(); ++i_node) {
+ for (size_t i_node_face = node_to_duplicate_face_row[i_node];
+ i_node_face < node_to_duplicate_face_row[i_node + 1] - 1; ++i_node_face) {
+ for (size_t j_node_face = i_node_face + 1; j_node_face < node_to_duplicate_face_row[i_node + 1]; ++j_node_face) {
+ if (dup_face_to_face[node_duplicate_face_list[i_node_face]] ==
+ dup_face_to_face[node_duplicate_face_list[j_node_face]])
+ continue;
+ if (is_same_face(node_duplicate_face_list[i_node_face], node_duplicate_face_list[j_node_face])) {
+ dup_face_to_face[node_duplicate_face_list[j_node_face]] =
+ dup_face_to_face[node_duplicate_face_list[i_node_face]];
+ nb_dup_faces++;
+ }
+ }
}
}
+ // compute face_id
+ Array<FaceId> dup_face_to_face_id(total_number_of_faces);
{
- descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
- for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
- descriptor.cell_face_is_reversed_vector[j] = Array<bool>(cell_nb_faces[j]);
+ FaceId face_id = 0;
+ for (size_t i_dup_face = 0; i_dup_face < total_number_of_faces; ++i_dup_face) {
+ if (dup_face_to_face[i_dup_face] == i_dup_face) {
+ dup_face_to_face_id[i_dup_face] = face_id;
+ ++face_id;
+ } else {
+ size_t i_face = dup_face_to_face[i_dup_face];
+ while (i_face != dup_face_to_face[i_face]) {
+ i_face = dup_face_to_face[i_face];
+ }
+ dup_face_to_face_id[i_dup_face] = dup_face_to_face_id[i_face];
+ }
}
- 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];
- descriptor.cell_face_is_reversed_vector[cell_number][cell_local_face] = reversed;
+ }
+
+ Array<unsigned int> node_to_face_row(node_to_duplicate_face_row.size());
+ {
+ unsigned int nb_faces = 0;
+ for (size_t i_node = 0; i_node < node_to_duplicate_face_row.size() - 1; ++i_node) {
+ node_to_face_row[i_node] = nb_faces;
+ for (size_t i_face = node_to_duplicate_face_row[i_node]; i_face < node_to_duplicate_face_row[i_node + 1];
+ ++i_face) {
+ if (dup_face_to_face[node_duplicate_face_list[i_face]] == node_duplicate_face_list[i_face]) {
+ ++nb_faces;
+ }
+ }
+ }
+ node_to_face_row[node_to_duplicate_face_row.size() - 1] = nb_faces;
+ }
+
+ Array<unsigned int> node_to_face_list(node_to_face_row[node_to_duplicate_face_row.size() - 1]);
+ {
+ unsigned int i_node_to_face = 0;
+ for (size_t i_node = 0; i_node < node_to_duplicate_face_row.size() - 1; ++i_node) {
+ for (size_t i_face = node_to_duplicate_face_row[i_node]; i_face < node_to_duplicate_face_row[i_node + 1];
+ ++i_face) {
+ if (dup_face_to_face[node_duplicate_face_list[i_face]] == node_duplicate_face_list[i_face]) {
+ node_to_face_list[i_node_to_face++] = dup_face_to_face_id[node_duplicate_face_list[i_face]];
+ }
+ }
+ }
+ }
+
+ descriptor.setNodeToFaceMatrix(ConnectivityMatrix(node_to_face_row, node_to_face_list));
+
+ Array<unsigned int> cell_to_face_row(cell_nb_faces.size() + 1);
+ cell_to_face_row[0] = 0;
+ for (size_t i = 0; i < cell_nb_faces.size(); ++i) {
+ cell_to_face_row[i + 1] = cell_to_face_row[i] + cell_nb_faces[i];
+ }
+
+ Array<unsigned int> cell_to_face_list(cell_to_face_row[cell_to_face_row.size() - 1]);
+ {
+ size_t i_cell_face = 0;
+ for (CellId cell_id = 0; cell_id < cell_nb_faces.size(); ++cell_id) {
+ for (size_t i_face = 0; i_face < cell_nb_faces[cell_id]; ++i_face) {
+ cell_to_face_list[i_cell_face++] = dup_face_to_face_id[cell_to_face_row[cell_id] + i_face];
}
}
}
+ descriptor.setCellToFaceMatrix(ConnectivityMatrix(cell_to_face_row, cell_to_face_list));
+
+ descriptor.setFaceNumberVector([&] {
+ Array<int> face_number_vector(total_number_of_faces - nb_dup_faces);
+ parallel_for(
+ face_number_vector.size(), PUGS_LAMBDA(const size_t l) { face_number_vector[l] = l; });
+ return face_number_vector;
+ }());
+
+ Array<unsigned int> face_ending(total_number_of_faces - nb_dup_faces + 1);
{
- 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.nodeIdList();
- ++l;
+ size_t i_face = 0;
+ face_ending[0] = 0;
+ for (size_t i_dup_face = 0; i_dup_face < dup_face_to_face.size(); ++i_dup_face) {
+ if (dup_face_to_face[i_dup_face] == i_dup_face) {
+ face_ending[i_face + 1] =
+ dup_face_to_node_row[i_dup_face + 1] - dup_face_to_node_row[i_dup_face] + face_ending[i_face];
+ ++i_face;
+ }
}
}
+ Array<unsigned int> faces_node_list(face_ending[face_ending.size() - 1]);
{
- // Face numbers may change if numbers are provided in the file
- 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;
+ size_t i_face = 0;
+ face_ending[0] = 0;
+ for (size_t i_dup_face = 0; i_dup_face < dup_face_to_face.size(); ++i_dup_face) {
+ if (dup_face_to_face[i_dup_face] == i_dup_face) {
+ size_t dup_face_node_begin = dup_face_to_node_row[i_dup_face];
+ size_t face_node_begin = face_ending[i_face];
+
+ for (size_t i_face_node = 0;
+ i_face_node < dup_face_to_node_row[i_dup_face + 1] - dup_face_to_node_row[i_dup_face]; ++i_face_node) {
+ faces_node_list[face_node_begin + i_face_node] = dup_faces_to_node_list[dup_face_node_begin + i_face_node];
+ }
+ ++i_face;
+ }
}
}
+
+ descriptor.setFaceToNodeMatrix(ConnectivityMatrix(face_ending, faces_node_list));
+
+ descriptor.setCellFaceIsReversed(cell_face_is_reversed);
}
template <size_t Dimension>
@@ -257,206 +598,331 @@ void
ConnectivityBuilderBase::_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()));
+
+ const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
+
+ Array<const unsigned int> face_to_edge_row = face_to_node_matrix.rowsMap();
+
+ const size_t total_number_of_face_edges = face_to_edge_row[face_to_edge_row.size() - 1];
+
+ const size_t total_number_of_node_by_face_edges = 2 * total_number_of_face_edges;
+
+ Array<unsigned int> face_edge_to_node_list(total_number_of_node_by_face_edges);
+ {
+ size_t i_edge_node = 0;
+ for (size_t i_face = 0; i_face < face_to_edge_row.size() - 1; ++i_face) {
+ const auto& face_node_list = face_to_node_matrix[i_face];
+ for (size_t i_node = 0; i_node < face_node_list.size() - 1; ++i_node) {
+ face_edge_to_node_list[i_edge_node++] = face_node_list[i_node];
+ face_edge_to_node_list[i_edge_node++] = face_node_list[i_node + 1];
+ }
+ face_edge_to_node_list[i_edge_node++] = face_node_list[face_node_list.size() - 1];
+ face_edge_to_node_list[i_edge_node++] = face_node_list[0];
}
}
- 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]);
+ auto node_number_vector = descriptor.nodeNumberVector();
+ Array<bool> face_edge_is_reversed(total_number_of_face_edges);
+ for (size_t i_edge = 0; i_edge < total_number_of_face_edges; ++i_edge) {
+ if (node_number_vector[face_edge_to_node_list[2 * i_edge]] >
+ node_number_vector[face_edge_to_node_list[2 * i_edge + 1]]) {
+ std::swap(face_edge_to_node_list[2 * i_edge], face_edge_to_node_list[2 * i_edge + 1]);
+ face_edge_is_reversed[i_edge] = true;
+ } else {
+ face_edge_is_reversed[i_edge] = false;
}
- 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;
+ }
+
+ const size_t total_number_of_edges = face_edge_is_reversed.size();
+ Array<unsigned int> node_to_face_edge_row(node_number_vector.size() + 1);
+ node_to_face_edge_row.fill(0);
+ for (size_t i_edge = 0; i_edge < total_number_of_edges; ++i_edge) {
+ for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
+ node_to_face_edge_row[face_edge_to_node_list[2 * i_edge + i_edge_node] + 1]++;
}
}
+ for (size_t i_node = 1; i_node < node_to_face_edge_row.size(); ++i_node) {
+ node_to_face_edge_row[i_node] += node_to_face_edge_row[i_node - 1];
+ }
+
+ Array<unsigned int> node_to_face_edge_list(node_to_face_edge_row[node_to_face_edge_row.size() - 1]);
{
- 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]);
+ Array<unsigned int> node_to_face_edge_row_idx(node_number_vector.size());
+ node_to_face_edge_row_idx.fill(0);
+
+ for (size_t i_edge = 0; i_edge < total_number_of_edges; ++i_edge) {
+ for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
+ const size_t node_id = face_edge_to_node_list[2 * i_edge + i_edge_node];
+
+ node_to_face_edge_list[node_to_face_edge_row[node_id] + node_to_face_edge_row_idx[node_id]] = i_edge;
+ node_to_face_edge_row_idx[node_id]++;
+ }
}
- 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;
+ }
+
+ Array<unsigned int> face_edge_to_edge(total_number_of_edges);
+ parallel_for(
+ total_number_of_edges, PUGS_LAMBDA(size_t i_edge) { face_edge_to_edge[i_edge] = i_edge; });
+
+ auto is_same_face = [=](const size_t i_edge, const size_t j_edge) {
+ auto i_edge_first_node = 2 * i_edge;
+ auto j_edge_first_node = 2 * j_edge;
+ return ((face_edge_to_node_list[i_edge_first_node] == face_edge_to_node_list[j_edge_first_node]) and
+ (face_edge_to_node_list[i_edge_first_node + 1] == face_edge_to_node_list[j_edge_first_node + 1]));
+ };
+
+ size_t nb_duplicate_edges = 0;
+ for (size_t i_node = 0; i_node < node_number_vector.size(); ++i_node) {
+ for (size_t i_node_face = node_to_face_edge_row[i_node]; i_node_face < node_to_face_edge_row[i_node + 1] - 1;
+ ++i_node_face) {
+ const unsigned int i_edge = node_to_face_edge_list[i_node_face];
+ const unsigned int i_edge_id = face_edge_to_edge[i_edge];
+ for (size_t j_node_face = i_node_face + 1; j_node_face < node_to_face_edge_row[i_node + 1]; ++j_node_face) {
+ const unsigned int j_edge = node_to_face_edge_list[j_node_face];
+ unsigned int& j_edge_id = face_edge_to_edge[j_edge];
+ if (i_edge_id != j_edge_id) {
+ if (is_same_face(i_edge, j_edge)) {
+ j_edge_id = i_edge_id;
+ nb_duplicate_edges++;
+ }
+ }
}
}
}
+ // compute edge_id
+ Array<EdgeId> dup_edge_to_edge_id(total_number_of_edges);
{
- 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;
+ EdgeId edge_id = 0;
+ for (size_t i_dup_edge = 0; i_dup_edge < total_number_of_edges; ++i_dup_edge) {
+ if (face_edge_to_edge[i_dup_edge] == i_dup_edge) {
+ dup_edge_to_edge_id[i_dup_edge] = edge_id++;
+ } else {
+ size_t i_edge = i_dup_edge;
+ do {
+ i_edge = face_edge_to_edge[i_edge];
+ } while (i_edge != face_edge_to_edge[i_edge]);
+
+ dup_edge_to_edge_id[i_dup_edge] = dup_edge_to_edge_id[i_edge];
+ }
}
}
+ Array<unsigned int> edge_to_node_list(2 * (total_number_of_edges - nb_duplicate_edges));
{
- // 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;
+ for (size_t i_dup_edge = 0; i_dup_edge < total_number_of_edges; ++i_dup_edge) {
+ if (face_edge_to_edge[i_dup_edge] == i_dup_edge) {
+ const EdgeId edge_id = dup_edge_to_edge_id[i_dup_edge];
+ edge_to_node_list[2 * edge_id] = face_edge_to_node_list[2 * i_dup_edge];
+ edge_to_node_list[2 * edge_id + 1] = face_edge_to_node_list[2 * i_dup_edge + 1];
+ }
}
}
- {
- descriptor.cell_to_edge_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];
+ descriptor.setEdgeNumberVector([&] {
+ Array<int> edge_number_vector(total_number_of_edges - nb_duplicate_edges);
+ parallel_for(
+ edge_number_vector.size(), PUGS_LAMBDA(const size_t i_edge) { edge_number_vector[i_edge] = i_edge; });
+ return edge_number_vector;
+ }());
- 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()) {
- throw NormalError("could not find this edge");
+ Array<unsigned int> edge_to_node_row(total_number_of_edges - nb_duplicate_edges + 1);
+ parallel_for(
+ edge_to_node_row.size(), PUGS_LAMBDA(const size_t i_edge) { edge_to_node_row[i_edge] = 2 * i_edge; });
+
+ descriptor.setEdgeToNodeMatrix(ConnectivityMatrix(edge_to_node_row, edge_to_node_list));
+
+ // Use real edge ids
+ for (size_t i_edge = 0; i_edge < face_edge_to_edge.size(); ++i_edge) {
+ face_edge_to_edge[i_edge] = dup_edge_to_edge_id[face_edge_to_edge[i_edge]];
+ }
+
+ descriptor.setFaceToEdgeMatrix(ConnectivityMatrix(face_to_edge_row, face_edge_to_edge));
+ descriptor.setFaceEdgeIsReversed(face_edge_is_reversed);
+
+ Array<size_t> node_to_duplicated_edge_id_list(node_to_face_edge_list.size());
+ for (size_t i_node_edge = 0; i_node_edge < node_to_duplicated_edge_id_list.size(); ++i_node_edge) {
+ node_to_duplicated_edge_id_list[i_node_edge] = dup_edge_to_edge_id[node_to_face_edge_list[i_node_edge]];
+ }
+
+ Array<bool> node_to_edge_is_duplicated(node_to_face_edge_list.size());
+ node_to_edge_is_duplicated.fill(false);
+
+ Array<unsigned int> node_nb_edges(node_number_vector.size());
+ node_nb_edges.fill(0);
+
+ for (size_t node_id = 0; node_id < node_number_vector.size(); ++node_id) {
+ size_t nb_dup_edges = 0;
+ for (EdgeId i_edge = node_to_face_edge_row[node_id]; i_edge < node_to_face_edge_row[node_id + 1] - 1; ++i_edge) {
+ if (not node_to_edge_is_duplicated[i_edge]) {
+ for (EdgeId j_edge = i_edge + 1; j_edge < node_to_face_edge_row[node_id + 1]; ++j_edge) {
+ if (node_to_duplicated_edge_id_list[i_edge] == node_to_duplicated_edge_id_list[j_edge]) {
+ node_to_edge_is_duplicated[j_edge] = true;
+ nb_dup_edges++;
}
- 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()) {
- throw NormalError("could not find this edge");
- }
- cell_edge_vector.push_back(i->second);
- }
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
- break;
+ }
+ node_nb_edges[node_id] = node_to_face_edge_row[node_id + 1] - node_to_face_edge_row[node_id] - nb_dup_edges;
+ }
+
+ Array<unsigned int> node_to_edge_row(node_number_vector.size() + 1);
+ node_to_edge_row[0] = 0;
+ for (size_t node_id = 0; node_id < node_number_vector.size(); ++node_id) {
+ node_to_edge_row[node_id + 1] = node_to_edge_row[node_id] + node_nb_edges[node_id];
+ }
+
+ Array<unsigned int> node_to_edge_list(node_to_edge_row[node_to_edge_row.size() - 1]);
+ {
+ size_t l = 0;
+ for (size_t i_edge_id = 0; i_edge_id < node_to_duplicated_edge_id_list.size(); ++i_edge_id) {
+ if (not node_to_edge_is_duplicated[i_edge_id]) {
+ node_to_edge_list[l++] = node_to_duplicated_edge_id_list[i_edge_id];
}
- case CellType::Prism: {
- constexpr int local_edge[12][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
- std::vector<unsigned int> cell_edge_vector;
- cell_edge_vector.reserve(9);
- for (int i_edge = 0; i_edge < 9; ++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()) {
- throw NormalError("could not find this edge");
- }
- cell_edge_vector.push_back(i->second);
+ }
+ }
+
+ descriptor.setNodeToEdgeMatrix(ConnectivityMatrix(node_to_edge_row, node_to_edge_list));
+
+ auto find_edge = [=](const NodeId& node0_id, const NodeId& node1_id) -> EdgeId {
+ auto [first_node_id, second_node_id] = [&] {
+ if (node_number_vector[node0_id] < node_number_vector[node1_id]) {
+ return std::make_pair(node0_id, node1_id);
+ } else {
+ return std::make_pair(node1_id, node0_id);
+ }
+ }();
+
+ for (size_t i_node_edge = node_to_edge_row[first_node_id]; i_node_edge < node_to_edge_row[first_node_id + 1];
+ ++i_node_edge) {
+ EdgeId edge_id = node_to_edge_list[i_node_edge];
+ if (edge_to_node_list[2 * edge_id + 1] == second_node_id) {
+ return edge_id;
+ }
+ }
+ throw UnexpectedError("Cannot find cell edge in edge list");
+ };
+
+ const auto& cell_to_node_matrix = descriptor.cellToNodeMatrix();
+ const auto& cell_type_vector = descriptor.cellTypeVector();
+ {
+ Array<unsigned int> cell_to_edge_row(cell_to_node_matrix.numberOfRows() + 1);
+ {
+ cell_to_edge_row[0] = 0;
+
+ for (CellId cell_id = 0; cell_id < cell_to_node_matrix.numberOfRows(); ++cell_id) {
+ switch (cell_type_vector[cell_id]) {
+ case CellType::Tetrahedron: {
+ cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 6;
+ break;
+ }
+ case CellType::Hexahedron: {
+ cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 12;
+ break;
+ }
+ case CellType::Prism: {
+ cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 9;
+ break;
+ }
+ case CellType::Pyramid: {
+ cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 2 * (cell_to_node_matrix[cell_id].size() - 1);
+ break;
+ }
+ case CellType::Diamond: {
+ cell_to_edge_row[cell_id + 1] = cell_to_edge_row[cell_id] + 3 * (cell_to_node_matrix[cell_id].size() - 2);
+ break;
+ }
+ default: {
+ std::stringstream error_msg;
+ error_msg << name(cell_type_vector[cell_id]) << ": unexpected cell type in dimension 3";
+ throw NotImplementedError(error_msg.str());
+ }
}
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
- break;
}
- case CellType::Pyramid: {
- const size_t number_of_edges = 2 * cell_nodes.size();
- std::vector<unsigned int> base_nodes;
- std::copy_n(cell_nodes.begin(), cell_nodes.size() - 1, std::back_inserter(base_nodes));
-
- std::vector<unsigned int> cell_edge_vector;
- cell_edge_vector.reserve(number_of_edges);
- for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
- Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
+ }
+
+ Array<unsigned int> cell_to_edge_list(cell_to_edge_row[cell_to_edge_row.size() - 1]);
+ {
+ size_t i_cell_edge = 0;
+ for (CellId cell_id = 0; cell_id < cell_to_node_matrix.numberOfRows(); ++cell_id) {
+ const auto& cell_nodes = cell_to_node_matrix[cell_id];
+
+ switch (cell_type_vector[cell_id]) {
+ case CellType::Tetrahedron: {
+ constexpr int local_edge[6][2] = {{0, 1}, {0, 2}, {0, 3}, {1, 2}, {2, 3}, {3, 1}};
+ for (int i_edge = 0; i_edge < 6; ++i_edge) {
+ const auto& e = local_edge[i_edge];
+ cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
}
- cell_edge_vector.push_back(i->second);
+ 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}};
- const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
- for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
- Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
+ for (int i_edge = 0; i_edge < 12; ++i_edge) {
+ const auto& e = local_edge[i_edge];
+ cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
}
- cell_edge_vector.push_back(i->second);
+ break;
}
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
- break;
- }
- case CellType::Diamond: {
- const size_t number_of_edges = 3 * cell_nodes.size();
- std::vector<unsigned int> base_nodes;
- std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
-
- std::vector<unsigned int> cell_edge_vector;
- cell_edge_vector.reserve(number_of_edges);
- for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
- Edge edge{{base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
+ case CellType::Prism: {
+ constexpr int local_edge[9][2] = {{0, 1}, {1, 2}, {2, 0}, {3, 4}, {4, 5}, {5, 3}, {0, 3}, {1, 4}, {2, 5}};
+ for (int i_edge = 0; i_edge < 9; ++i_edge) {
+ const auto& e = local_edge[i_edge];
+ cell_to_edge_list[i_cell_edge++] = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
}
- cell_edge_vector.push_back(i->second);
+ break;
}
+ case CellType::Pyramid: {
+ auto base_nodes = [&](size_t i) { return cell_nodes[i]; };
- const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
- for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
- Edge edge{{base_nodes[i_edge], top_vertex}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
+ for (size_t i_edge = 0; i_edge < cell_nodes.size() - 1; ++i_edge) {
+ cell_to_edge_list[i_cell_edge++] =
+ find_edge(base_nodes(i_edge), base_nodes((i_edge + 1) % (cell_nodes.size() - 1)));
}
- cell_edge_vector.push_back(i->second);
- }
- const unsigned int bottom_vertex = cell_nodes[0];
- for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
- Edge edge{{base_nodes[i_edge], bottom_vertex}, node_number_vector};
- auto i = edge_id_map.find(edge);
- if (i == edge_id_map.end()) {
- throw NormalError("could not find this edge");
+ const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_edge = 0; i_edge < cell_nodes.size() - 1; ++i_edge) {
+ cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), top_vertex);
}
- cell_edge_vector.push_back(i->second);
+ break;
}
+ case CellType::Diamond: {
+ auto base_nodes = [&](size_t i) { return cell_nodes[i + 1]; };
- descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
+ for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+ cell_to_edge_list[i_cell_edge++] =
+ find_edge(base_nodes(i_edge), base_nodes((i_edge + 1) % (cell_nodes.size() - 2)));
+ }
- break;
- }
- default: {
- std::stringstream error_msg;
- error_msg << name(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 3";
- throw NotImplementedError(error_msg.str());
- }
+ {
+ const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+ cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), top_vertex);
+ }
+ }
+
+ {
+ const unsigned int bottom_vertex = cell_nodes[0];
+ for (size_t i_edge = 0; i_edge < cell_nodes.size() - 2; ++i_edge) {
+ cell_to_edge_list[i_cell_edge++] = find_edge(base_nodes(i_edge), bottom_vertex);
+ }
+ }
+
+ break;
+ }
+ default: {
+ std::stringstream error_msg;
+ error_msg << name(cell_type_vector[cell_id]) << ": unexpected cell type in dimension 3";
+ throw NotImplementedError(error_msg.str());
+ }
+ }
}
}
+
+ descriptor.setCellToEdgeMatrix(ConnectivityMatrix(cell_to_edge_row, cell_to_edge_list));
}
}
diff --git a/src/mesh/ConnectivityBuilderBase.hpp b/src/mesh/ConnectivityBuilderBase.hpp
index f9b0d28a7a9389093348648573c85875f1e4a7de..a621e91f46ffe7e6ecb5ea6d496934a03571ba6e 100644
--- a/src/mesh/ConnectivityBuilderBase.hpp
+++ b/src/mesh/ConnectivityBuilderBase.hpp
@@ -15,9 +15,6 @@ class ConnectivityDescriptor;
class ConnectivityBuilderBase
{
protected:
- template <size_t Dimension>
- class ConnectivityFace;
-
std::shared_ptr<const IConnectivity> m_connectivity;
template <size_t Dimension>
@@ -37,190 +34,4 @@ class ConnectivityBuilderBase
~ConnectivityBuilderBase() = default;
};
-template <>
-class ConnectivityBuilderBase::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;
- }
- };
-
- private:
- const std::vector<int>& m_node_number_vector;
-
- unsigned int m_node0_id;
- unsigned int m_node1_id;
-
- bool m_reversed;
-
- public:
- std::vector<unsigned int>
- nodeIdList() const
- {
- return {m_node0_id, m_node1_id};
- }
-
- bool
- reversed() const
- {
- return m_reversed;
- }
-
- PUGS_INLINE
- bool
- operator==(const ConnectivityFace& f) const
- {
- return ((m_node0_id == f.m_node0_id) and (m_node1_id == f.m_node1_id));
- }
-
- PUGS_INLINE
- bool
- operator<(const ConnectivityFace& f) const
- {
- 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])));
- }
-
- PUGS_INLINE
- ConnectivityFace(const std::vector<unsigned int>& node_id_list, const std::vector<int>& node_number_vector)
- : m_node_number_vector(node_number_vector)
- {
- 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;
- }
- }
-
- PUGS_INLINE
- ConnectivityFace(const ConnectivityFace&) = default;
-
- PUGS_INLINE
- ~ConnectivityFace() = default;
-};
-
-template <>
-class ConnectivityBuilderBase::ConnectivityFace<3>
-{
- public:
- 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;
- }
- };
-
- private:
- bool m_reversed;
- std::vector<NodeId::base_type> m_node_id_list;
- const std::vector<int>& m_node_number_vector;
-
- PUGS_INLINE
- 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;
- }
-
- public:
- PUGS_INLINE
- const bool&
- reversed() const
- {
- return m_reversed;
- }
-
- PUGS_INLINE
- const std::vector<unsigned int>&
- nodeIdList() const
- {
- return m_node_id_list;
- }
-
- PUGS_INLINE
- 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_number_vector(node_number_vector)
- {
- ;
- }
-
- 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;
- }
-
- PUGS_INLINE
- 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();
- }
-
- PUGS_INLINE
- ConnectivityFace(const ConnectivityFace&) = default;
-
- PUGS_INLINE
- ConnectivityFace() = delete;
-
- PUGS_INLINE
- ~ConnectivityFace() = default;
-};
-
#endif // CONNECTIVITY_BUILDER_BASE_HPP
diff --git a/src/mesh/ConnectivityComputer.cpp b/src/mesh/ConnectivityComputer.cpp
index aa8ca779cfb965ffb3e85e4feb1b4b9feabb79ca..1d5514627a5dac75d41ba5e4d08fa182d37e552f 100644
--- a/src/mesh/ConnectivityComputer.cpp
+++ b/src/mesh/ConnectivityComputer.cpp
@@ -7,62 +7,90 @@
template <typename ConnectivityType>
PUGS_INLINE ConnectivityMatrix
-ConnectivityComputer::computeConnectivityMatrix(const ConnectivityType& connectivity,
- ItemType item_type,
- ItemType child_item_type) const
+ConnectivityComputer::computeInverseConnectivityMatrix(const ConnectivityType& connectivity,
+ ItemType item_type,
+ 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);
+ if (item_type < child_item_type) {
+ 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);
+
+ switch (child_item_type) {
+ case ItemType::cell: {
+ item_to_child_item_matrix =
+ this->_computeInverseConnectivity<ConnectivityType, ItemType::cell>(connectivity, child_to_item_matrix);
+ break;
+ }
+ case ItemType::face: {
+ item_to_child_item_matrix =
+ this->_computeInverseConnectivity<ConnectivityType, ItemType::face>(connectivity, child_to_item_matrix);
+ break;
+ }
+ case ItemType::edge: {
+ item_to_child_item_matrix =
+ this->_computeInverseConnectivity<ConnectivityType, ItemType::edge>(connectivity, child_to_item_matrix);
+ break;
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("node cannot be child item when computing inverse connectivity");
+ }
+ // LCOV_EXCL_STOP
+ }
+ } else {
+ std::stringstream error_msg;
+ error_msg << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type);
+ throw UnexpectedError(error_msg.str());
+ }
- item_to_child_item_matrix = this->_computeInverse(child_to_item_matrix);
+ return item_to_child_item_matrix;
} else {
std::stringstream error_msg;
- error_msg << "unable to compute connectivity " << itemName(item_type) << " -> " << itemName(child_item_type);
+ error_msg << "cannot deduce " << itemName(item_type) << " -> " << itemName(child_item_type) << " connectivity";
throw UnexpectedError(error_msg.str());
}
-
- return item_to_child_item_matrix;
}
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity1D&,
- ItemType,
- ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity1D&,
+ ItemType,
+ ItemType) const;
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity2D&,
- ItemType,
- ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity2D&,
+ ItemType,
+ ItemType) const;
-template ConnectivityMatrix ConnectivityComputer::computeConnectivityMatrix(const Connectivity3D&,
- ItemType,
- ItemType) const;
+template ConnectivityMatrix ConnectivityComputer::computeInverseConnectivityMatrix(const Connectivity3D&,
+ ItemType,
+ ItemType) const;
+template <typename ConnectivityType, ItemType child_item_type>
ConnectivityMatrix
-ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_matrix) const
+ConnectivityComputer::_computeInverseConnectivity(const ConnectivityType& connectivity,
+ const ConnectivityMatrix& child_to_item_matrix) const
{
- const size_t& number_of_rows = item_to_child_matrix.numberOfRows();
+ const size_t number_of_transposed_columns = child_to_item_matrix.numberOfRows();
- if ((item_to_child_matrix.values().size() > 0)) {
- const size_t& number_of_columns = max(item_to_child_matrix.values());
+ if ((child_to_item_matrix.values().size() > 0)) {
+ const size_t number_of_transposed_rows = max(child_to_item_matrix.values()) + 1;
- Array<uint32_t> transposed_next_free_column_index(number_of_columns + 1);
+ Array<uint32_t> transposed_next_free_column_index(number_of_transposed_rows);
transposed_next_free_column_index.fill(0);
Array<uint32_t> transposed_rows_map(transposed_next_free_column_index.size() + 1);
transposed_rows_map.fill(0);
- for (size_t i = 0; i < number_of_rows; ++i) {
- for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
- transposed_rows_map[item_to_child_matrix.values()[j] + 1]++;
+ for (size_t i = 0; i < number_of_transposed_columns; ++i) {
+ for (size_t j = child_to_item_matrix.rowsMap()[i]; j < child_to_item_matrix.rowsMap()[i + 1]; ++j) {
+ transposed_rows_map[child_to_item_matrix.values()[j] + 1]++;
}
}
for (size_t i = 1; i < transposed_rows_map.size(); ++i) {
transposed_rows_map[i] += transposed_rows_map[i - 1];
}
Array<uint32_t> transposed_column_indices(transposed_rows_map[transposed_rows_map.size() - 1]);
-
- for (size_t i = 0; i < number_of_rows; ++i) {
- for (size_t j = item_to_child_matrix.rowsMap()[i]; j < item_to_child_matrix.rowsMap()[i + 1]; ++j) {
- size_t i_column_index = item_to_child_matrix.values()[j];
+ for (size_t i = 0; i < number_of_transposed_columns; ++i) {
+ for (size_t j = child_to_item_matrix.rowsMap()[i]; j < child_to_item_matrix.rowsMap()[i + 1]; ++j) {
+ size_t i_column_index = child_to_item_matrix.values()[j];
uint32_t& shift = transposed_next_free_column_index[i_column_index];
transposed_column_indices[transposed_rows_map[i_column_index] + shift] = i;
@@ -71,6 +99,18 @@ ConnectivityComputer::_computeInverse(const ConnectivityMatrix& item_to_child_ma
}
}
+ auto target_item_number = connectivity.template number<child_item_type>();
+ // Finally one sorts target item_ids for parallel reproducibility
+ for (size_t i = 0; i < number_of_transposed_rows; ++i) {
+ auto row_begining = &(transposed_column_indices[transposed_rows_map[i]]);
+ auto row_size = (transposed_rows_map[i + 1] - transposed_rows_map[i]);
+ std::sort(row_begining, row_begining + row_size,
+ [&target_item_number](const ItemIdT<child_item_type> item0_id,
+ const ItemIdT<child_item_type> item1_id) {
+ return target_item_number[item0_id] < target_item_number[item1_id];
+ });
+ }
+
return ConnectivityMatrix{transposed_rows_map, transposed_column_indices};
} else {
// empty connectivity
diff --git a/src/mesh/ConnectivityComputer.hpp b/src/mesh/ConnectivityComputer.hpp
index b033a25d5ccc04db84c194e79e492a081216a49e..7093aeef17707c37452b5e6874f6ed50a285dbb3 100644
--- a/src/mesh/ConnectivityComputer.hpp
+++ b/src/mesh/ConnectivityComputer.hpp
@@ -7,13 +7,15 @@
class ConnectivityComputer
{
private:
- ConnectivityMatrix _computeInverse(const ConnectivityMatrix& item_to_child_matrix) const;
+ template <typename ConnectivityType, ItemType child_item_type>
+ ConnectivityMatrix _computeInverseConnectivity(const ConnectivityType& connectivity,
+ const ConnectivityMatrix& item_to_child_matrix) const;
public:
template <typename ConnectivityType>
- ConnectivityMatrix computeConnectivityMatrix(const ConnectivityType& connectivity,
- ItemType item_type,
- ItemType child_item_type) const;
+ ConnectivityMatrix computeInverseConnectivityMatrix(const ConnectivityType& connectivity,
+ ItemType item_type,
+ ItemType child_item_type) const;
template <typename ItemOfItem, typename ConnectivityType>
void computeLocalItemNumberInChildItem(const ConnectivityType& connectivity) const;
diff --git a/src/mesh/ConnectivityDescriptor.hpp b/src/mesh/ConnectivityDescriptor.hpp
index 8f26988eadaa01dfec8aa0f14404a82819273475..dacf84d3cc84ab13629c09f163a16739db70dba3 100644
--- a/src/mesh/ConnectivityDescriptor.hpp
+++ b/src/mesh/ConnectivityDescriptor.hpp
@@ -2,6 +2,7 @@
#define CONNECTIVITY_DESCRIPTOR_HPP
#include <mesh/CellType.hpp>
+#include <mesh/ConnectivityMatrix.hpp>
#include <mesh/ItemOfItemType.hpp>
#include <mesh/RefItemList.hpp>
#include <utils/PugsTraits.hpp>
@@ -16,68 +17,337 @@ class ConnectivityDescriptor
std::vector<RefEdgeList> m_ref_edge_list_vector;
std::vector<RefNodeList> m_ref_node_list_vector;
+ ConnectivityMatrix m_cell_to_face_matrix = ConnectivityMatrix{true};
+ ConnectivityMatrix m_cell_to_edge_matrix = ConnectivityMatrix{true};
+ ConnectivityMatrix m_cell_to_node_matrix = ConnectivityMatrix{true};
+
+ ConnectivityMatrix m_face_to_edge_matrix = ConnectivityMatrix{true};
+ ConnectivityMatrix m_face_to_node_matrix = ConnectivityMatrix{true};
+
+ ConnectivityMatrix m_edge_to_node_matrix = ConnectivityMatrix{true};
+
+ ConnectivityMatrix m_node_to_face_matrix = ConnectivityMatrix{true};
+ ConnectivityMatrix m_node_to_edge_matrix = ConnectivityMatrix{true};
+
+ Array<const bool> m_cell_face_is_reversed;
+ Array<const bool> m_face_edge_is_reversed;
+
+ Array<const CellType> m_cell_type_vector;
+
+ Array<const int> m_cell_number_vector;
+ Array<const int> m_face_number_vector;
+ Array<const int> m_edge_number_vector;
+ Array<const int> m_node_number_vector;
+
+ Array<const int> m_cell_owner_vector;
+ Array<const int> m_face_owner_vector;
+ Array<const int> m_edge_owner_vector;
+ Array<const int> m_node_owner_vector;
+
public:
- std::vector<std::vector<unsigned int>> cell_to_node_vector;
- std::vector<std::vector<unsigned int>> cell_to_face_vector;
- std::vector<std::vector<unsigned int>> cell_to_edge_vector;
+ void
+ setCellNumberVector(const Array<const int>& cell_number_vector)
+ {
+ // No check since it can change reading file for instance
+ m_cell_number_vector = cell_number_vector;
+ }
- std::vector<std::vector<unsigned int>> face_to_node_vector;
- std::vector<std::vector<unsigned int>> face_to_edge_vector;
+ PUGS_INLINE
+ const Array<const int>&
+ cellNumberVector() const
+ {
+ return m_cell_number_vector;
+ }
- std::vector<std::vector<unsigned int>> edge_to_node_vector;
+ void
+ setFaceNumberVector(const Array<const int>& face_number_vector)
+ {
+ // No check since it can change reading file for instance
+ m_face_number_vector = face_number_vector;
+ }
- template <typename ItemOfItemT>
- auto&
- itemOfItemVector()
+ PUGS_INLINE
+ const Array<const int>&
+ faceNumberVector() const
{
- if constexpr (std::is_same_v<ItemOfItemT, NodeOfCell>) {
- 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");
- }
+ return m_face_number_vector;
+ }
+
+ void
+ setEdgeNumberVector(const Array<const int>& edge_number_vector)
+ {
+ // No check since it can change reading file for instance
+ m_edge_number_vector = edge_number_vector;
}
- std::vector<Array<bool>> cell_face_is_reversed_vector;
- std::vector<Array<bool>> face_edge_is_reversed_vector;
+ PUGS_INLINE
+ const Array<const int>&
+ edgeNumberVector() const
+ {
+ return m_edge_number_vector;
+ }
- std::vector<CellType> cell_type_vector;
+ void
+ setNodeNumberVector(const Array<const int>& node_number_vector)
+ {
+ // No check since it can change reading file for instance
+ m_node_number_vector = node_number_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;
+ PUGS_INLINE
+ const Array<const int>&
+ nodeNumberVector() const
+ {
+ return m_node_number_vector;
+ }
template <ItemType item_type>
- const std::vector<int>&
+ Array<const int>
itemNumberVector() const
{
if constexpr (item_type == ItemType::cell) {
- return cell_number_vector;
+ return m_cell_number_vector;
} else if constexpr (item_type == ItemType::face) {
- return face_number_vector;
+ return m_face_number_vector;
} else if constexpr (item_type == ItemType::edge) {
- return edge_number_vector;
+ return m_edge_number_vector;
} else if constexpr (item_type == ItemType::node) {
- return node_number_vector;
+ return m_node_number_vector;
} else {
static_assert(is_false_item_type_v<item_type>, "Unexpected item type");
}
}
- std::vector<int> cell_owner_vector;
- std::vector<int> face_owner_vector;
- std::vector<int> edge_owner_vector;
- std::vector<int> node_owner_vector;
+ void
+ setCellOwnerVector(const Array<const int>& cell_owner_vector)
+ {
+ Assert(m_cell_owner_vector.size() == 0);
+ m_cell_owner_vector = cell_owner_vector;
+ }
+
+ PUGS_INLINE
+ const Array<const int>&
+ cellOwnerVector() const
+ {
+ return m_cell_owner_vector;
+ }
+
+ void
+ setFaceOwnerVector(const Array<const int>& face_owner_vector)
+ {
+ Assert(m_face_owner_vector.size() == 0);
+ m_face_owner_vector = face_owner_vector;
+ }
+
+ PUGS_INLINE
+ const Array<const int>&
+ faceOwnerVector() const
+ {
+ return m_face_owner_vector;
+ }
+
+ void
+ setEdgeOwnerVector(const Array<const int>& edge_owner_vector)
+ {
+ Assert(m_edge_owner_vector.size() == 0);
+ m_edge_owner_vector = edge_owner_vector;
+ }
+
+ PUGS_INLINE
+ const Array<const int>&
+ edgeOwnerVector() const
+ {
+ return m_edge_owner_vector;
+ }
+
+ void
+ setNodeOwnerVector(const Array<const int>& node_owner_vector)
+ {
+ Assert(m_node_owner_vector.size() == 0);
+ m_node_owner_vector = node_owner_vector;
+ }
+
+ PUGS_INLINE
+ const Array<const int>&
+ nodeOwnerVector() const
+ {
+ return m_node_owner_vector;
+ }
+
+ void
+ setCellFaceIsReversed(const Array<const bool>& cell_face_is_reversed)
+ {
+ Assert(m_cell_face_is_reversed.size() == 0);
+ m_cell_face_is_reversed = cell_face_is_reversed;
+ }
+
+ PUGS_INLINE
+ const Array<const bool>&
+ cellFaceIsReversed() const
+ {
+ return m_cell_face_is_reversed;
+ }
+
+ void
+ setFaceEdgeIsReversed(const Array<const bool>& face_edge_is_reversed)
+ {
+ Assert(m_face_edge_is_reversed.size() == 0);
+ m_face_edge_is_reversed = face_edge_is_reversed;
+ }
+
+ PUGS_INLINE
+ const Array<const bool>&
+ faceEdgeIsReversed() const
+ {
+ return m_face_edge_is_reversed;
+ }
+
+ void
+ setCellTypeVector(const Array<const CellType>& cell_type_vector)
+ {
+ Assert(m_face_edge_is_reversed.size() == 0);
+ m_cell_type_vector = cell_type_vector;
+ }
+
+ PUGS_INLINE
+ const Array<const CellType>&
+ cellTypeVector() const
+ {
+ return m_cell_type_vector;
+ }
+
+ void
+ setCellToFaceMatrix(const ConnectivityMatrix& cell_to_face_matrix)
+ {
+ Assert(m_cell_to_face_matrix.numberOfRows() == 0);
+ m_cell_to_face_matrix = cell_to_face_matrix;
+ }
+
+ void
+ setCellToEdgeMatrix(const ConnectivityMatrix& cell_to_edge_matrix)
+ {
+ Assert(m_cell_to_edge_matrix.numberOfRows() == 0);
+ m_cell_to_edge_matrix = cell_to_edge_matrix;
+ }
+
+ void
+ setCellToNodeMatrix(const ConnectivityMatrix& cell_to_node_matrix)
+ {
+ Assert(m_cell_to_node_matrix.numberOfRows() == 0);
+ m_cell_to_node_matrix = cell_to_node_matrix;
+ }
+
+ void
+ setFaceToEdgeMatrix(const ConnectivityMatrix& face_to_edge_matrix)
+ {
+ Assert(m_face_to_edge_matrix.numberOfRows() == 0);
+ m_face_to_edge_matrix = face_to_edge_matrix;
+ }
+
+ void
+ setFaceToNodeMatrix(const ConnectivityMatrix& face_to_node_matrix)
+ {
+ Assert(m_face_to_node_matrix.numberOfRows() == 0);
+ m_face_to_node_matrix = face_to_node_matrix;
+ }
+
+ void
+ setEdgeToNodeMatrix(const ConnectivityMatrix& edge_to_node_matrix)
+ {
+ Assert(m_edge_to_node_matrix.numberOfRows() == 0);
+ m_edge_to_node_matrix = edge_to_node_matrix;
+ }
+
+ void
+ setNodeToFaceMatrix(const ConnectivityMatrix& node_to_face_matrix)
+ {
+ Assert(m_node_to_face_matrix.numberOfRows() == 0);
+ m_node_to_face_matrix = node_to_face_matrix;
+ }
+
+ void
+ setNodeToEdgeMatrix(const ConnectivityMatrix& node_to_edge_matrix)
+ {
+ Assert(m_node_to_edge_matrix.numberOfRows() == 0);
+ m_node_to_edge_matrix = node_to_edge_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ cellToFaceMatrix() const
+ {
+ return m_cell_to_face_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ cellToEdgeMatrix() const
+ {
+ return m_cell_to_edge_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ cellToNodeMatrix() const
+ {
+ return m_cell_to_node_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ faceToEdgeMatrix() const
+ {
+ return m_face_to_edge_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ faceToNodeMatrix() const
+ {
+ return m_face_to_node_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ edgeToNodeMatrix() const
+ {
+ return m_edge_to_node_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ nodeToFaceMatrix() const
+ {
+ return m_node_to_face_matrix;
+ }
+
+ PUGS_INLINE
+ const ConnectivityMatrix&
+ nodeToEdgeMatrix() const
+ {
+ return m_node_to_edge_matrix;
+ }
+
+ template <typename ItemOfItemT>
+ auto&
+ itemOfItemVector()
+ {
+ if constexpr (std::is_same_v<ItemOfItemT, NodeOfCell>) {
+ return m_cell_to_node_matrix;
+ } else if constexpr (std::is_same_v<ItemOfItemT, FaceOfCell>) {
+ return m_cell_to_face_matrix;
+ } else if constexpr (std::is_same_v<ItemOfItemT, EdgeOfCell>) {
+ return m_cell_to_edge_matrix;
+ } else if constexpr (std::is_same_v<ItemOfItemT, EdgeOfFace>) {
+ return m_face_to_edge_matrix;
+ } else if constexpr (std::is_same_v<ItemOfItemT, NodeOfFace>) {
+ return m_face_to_node_matrix;
+ } else if constexpr (std::is_same_v<ItemOfItemT, NodeOfEdge>) {
+ return m_edge_to_node_matrix;
+ } else {
+ static_assert(is_false_v<ItemOfItemT>, "Unexpected item of item type");
+ }
+ }
template <ItemType item_type>
const std::vector<RefItemList<item_type>>&
@@ -114,7 +384,7 @@ class ConnectivityDescriptor
}
ConnectivityDescriptor& operator=(const ConnectivityDescriptor&) = delete;
- ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
+ ConnectivityDescriptor& operator=(ConnectivityDescriptor&&) = delete;
ConnectivityDescriptor() = default;
ConnectivityDescriptor(const ConnectivityDescriptor&) = default;
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index fd5b58c6daaf5fd44c0cb440e7af472b9f9aebbf..641d7f8666946a448d98e0a3311d5bd2a94329fc 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -145,19 +145,19 @@ template <int Dimension>
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void
ConnectivityDispatcher<Dimension>::_gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
- std::vector<std::remove_const_t<DataType>>& gathered_vector)
+ Array<std::remove_const_t<DataType>>& gathered_array)
{
std::vector<Array<const DataType>> recv_item_data_by_proc = this->exchange(data_to_gather);
const auto& recv_id_correspondance_by_proc = this->_dispatchedInfo<item_type>().m_recv_id_correspondance_by_proc;
Assert(recv_id_correspondance_by_proc.size() == parallel::size());
- gathered_vector.resize(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
+ gathered_array = Array<std::remove_const_t<DataType>>(this->_dispatchedInfo<item_type>().m_number_to_id_map.size());
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
Assert(recv_id_correspondance_by_proc[i_rank].size() == recv_item_data_by_proc[i_rank].size());
for (size_t r = 0; r < recv_id_correspondance_by_proc[i_rank].size(); ++r) {
- const auto& item_id = recv_id_correspondance_by_proc[i_rank][r];
- gathered_vector[item_id] = recv_item_data_by_proc[i_rank][r];
+ const auto& item_id = recv_id_correspondance_by_proc[i_rank][r];
+ gathered_array[item_id] = recv_item_data_by_proc[i_rank][r];
}
}
}
@@ -167,7 +167,7 @@ template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
void
ConnectivityDispatcher<Dimension>::_gatherFrom(
const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
- std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector)
+ Array<std::remove_const_t<DataType>>& gathered_array)
{
using MutableDataType = std::remove_const_t<DataType>;
@@ -201,17 +201,23 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
parallel::exchange(data_to_send_by_proc, recv_data_to_gather_by_proc);
- const auto& item_list_to_recv_size_by_proc = this->_dispatchedInfo<item_type>().m_list_to_recv_size_by_proc;
+ const size_t recv_array_size = [&] {
+ size_t size = 0;
+ for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+ size += recv_data_to_gather_by_proc[i_rank].size();
+ }
+ return size;
+ }();
- for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
- int l = 0;
- for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
- Array<MutableDataType> data_vector(number_of_sub_item_per_item_to_recv_by_proc[i_rank][i]);
- for (size_t k = 0; k < data_vector.size(); ++k) {
- data_vector[k] = recv_data_to_gather_by_proc[i_rank][l++];
+ gathered_array = Array<std::remove_const_t<DataType>>(recv_array_size);
+ {
+ size_t l = 0;
+ for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
+ for (size_t j = 0; j < recv_data_to_gather_by_proc[i_rank].size(); ++j) {
+ gathered_array[l++] = recv_data_to_gather_by_proc[i_rank][j];
}
- gathered_vector.emplace_back(data_vector);
}
+ Assert(gathered_array.size() == l);
}
}
@@ -342,6 +348,8 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
const auto& recv_item_of_item_numbers_by_proc =
this->_dispatchedInfo<ItemOfItemT>().m_sub_item_numbers_to_recv_by_proc;
+ std::vector<std::vector<unsigned int>> item_to_subitem_legacy;
+ size_t number_of_node_by_cell = 0;
for (size_t i_rank = 0; i_rank < parallel::size(); ++i_rank) {
int l = 0;
for (size_t i = 0; i < item_list_to_recv_size_by_proc[i_rank]; ++i) {
@@ -351,9 +359,30 @@ ConnectivityDispatcher<Dimension>::_buildItemToSubItemDescriptor()
Assert(searched_sub_item_id != sub_item_number_id_map.end());
sub_item_vector.push_back(searched_sub_item_id->second);
}
- m_new_descriptor.itemOfItemVector<ItemOfItemT>().emplace_back(sub_item_vector);
+ number_of_node_by_cell += sub_item_vector.size();
+
+ item_to_subitem_legacy.emplace_back(sub_item_vector);
+ }
+ }
+ Array<unsigned int> item_to_subitem_row_map(item_to_subitem_legacy.size() + 1);
+ Array<unsigned int> item_to_subitem_list(number_of_node_by_cell);
+
+ item_to_subitem_row_map.fill(10000000);
+ item_to_subitem_list.fill(10000000);
+
+ item_to_subitem_row_map[0] = 0;
+ for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
+ item_to_subitem_row_map[i + 1] = item_to_subitem_row_map[i] + item_to_subitem_legacy[i].size();
+ }
+ size_t l = 0;
+ for (size_t i = 0; i < item_to_subitem_legacy.size(); ++i) {
+ const auto& subitem_list = item_to_subitem_legacy[i];
+ for (size_t j = 0; j < subitem_list.size(); ++j, ++l) {
+ item_to_subitem_list[l] = subitem_list[j];
}
}
+
+ m_new_descriptor.itemOfItemVector<ItemOfItemT>() = ConnectivityMatrix(item_to_subitem_row_map, item_to_subitem_list);
}
template <int Dimension>
@@ -582,7 +611,11 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
this->_buildSubItemNumberToIdMap<EdgeOfCell>();
this->_buildItemToExchangeLists<ItemType::edge>();
- this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), m_new_descriptor.edge_number_vector);
+ m_new_descriptor.setEdgeNumberVector([&] {
+ Array<int> edge_number_vector;
+ this->_gatherFrom(m_connectivity.template number<ItemType::edge>(), edge_number_vector);
+ return edge_number_vector;
+ }());
this->_buildItemToSubItemDescriptor<EdgeOfCell>();
@@ -594,9 +627,17 @@ ConnectivityDispatcher<Dimension>::_dispatchEdges()
this->_buildSubItemNumbersToRecvByProc<EdgeOfFace>();
this->_buildItemToSubItemDescriptor<EdgeOfFace>();
- this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), m_new_descriptor.face_edge_is_reversed_vector);
+ m_new_descriptor.setFaceEdgeIsReversed([&] {
+ Array<bool> face_edge_is_reversed;
+ this->_gatherFrom(m_connectivity.faceEdgeIsReversed(), face_edge_is_reversed);
+ return face_edge_is_reversed;
+ }());
- this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, m_new_descriptor.edge_owner_vector);
+ m_new_descriptor.setEdgeOwnerVector([&] {
+ Array<int> edge_owner_vector;
+ this->_gatherFrom(this->_dispatchedInfo<ItemType::edge>().m_new_owner, edge_owner_vector);
+ return edge_owner_vector;
+ }());
this->_buildItemReferenceList<ItemType::edge>();
}
@@ -616,13 +657,25 @@ ConnectivityDispatcher<Dimension>::_dispatchFaces()
this->_buildSubItemNumbersToRecvByProc<NodeOfFace>();
this->_buildItemToSubItemDescriptor<NodeOfFace>();
- this->_gatherFrom(m_connectivity.template number<ItemType::face>(), m_new_descriptor.face_number_vector);
+ m_new_descriptor.setFaceNumberVector([&] {
+ Array<int> face_number_vector;
+ this->_gatherFrom(m_connectivity.template number<ItemType::face>(), face_number_vector);
+ return face_number_vector;
+ }());
this->_buildItemToSubItemDescriptor<FaceOfCell>();
- this->_gatherFrom(m_connectivity.cellFaceIsReversed(), m_new_descriptor.cell_face_is_reversed_vector);
+ m_new_descriptor.setCellFaceIsReversed([&] {
+ Array<bool> cell_face_is_reversed;
+ this->_gatherFrom(m_connectivity.cellFaceIsReversed(), cell_face_is_reversed);
+ return cell_face_is_reversed;
+ }());
- this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, m_new_descriptor.face_owner_vector);
+ m_new_descriptor.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector;
+ this->_gatherFrom(this->_dispatchedInfo<ItemType::face>().m_new_owner, face_owner_vector);
+ return face_owner_vector;
+ }());
this->_buildItemReferenceList<ItemType::face>();
}
@@ -647,18 +700,39 @@ ConnectivityDispatcher<Dimension>::ConnectivityDispatcher(const ConnectivityType
this->_buildSubItemNumbersToRecvByProc<NodeOfCell>();
- this->_gatherFrom(m_connectivity.template number<ItemType::cell>(), m_new_descriptor.cell_number_vector);
+ m_new_descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector;
+ this->_gatherFrom(m_connectivity.template number<ItemType::cell>(), cell_number_vector);
+ return cell_number_vector;
+ }());
this->_buildSubItemNumberToIdMap<NodeOfCell>();
this->_buildItemToExchangeLists<ItemType::node>();
- // Fill new descriptor
- this->_gatherFrom(m_connectivity.cellType(), m_new_descriptor.cell_type_vector);
- this->_gatherFrom(this->_dispatchedInfo<ItemType::cell>().m_new_owner, m_new_descriptor.cell_owner_vector);
-
- this->_gatherFrom(m_connectivity.template number<ItemType::node>(), m_new_descriptor.node_number_vector);
- this->_gatherFrom(this->_dispatchedInfo<ItemType::node>().m_new_owner, m_new_descriptor.node_owner_vector);
+ m_new_descriptor.setCellTypeVector([&] {
+ Array<CellType> cell_type_vector;
+ this->_gatherFrom(m_connectivity.cellType(), cell_type_vector);
+ return cell_type_vector;
+ }());
+
+ m_new_descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector;
+ this->_gatherFrom(this->_dispatchedInfo<ItemType::cell>().m_new_owner, cell_owner_vector);
+ return cell_owner_vector;
+ }());
+
+ m_new_descriptor.setNodeNumberVector([&] {
+ Array<int> node_number_vector;
+ this->_gatherFrom(m_connectivity.template number<ItemType::node>(), node_number_vector);
+ return node_number_vector;
+ }());
+
+ m_new_descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector;
+ this->_gatherFrom(this->_dispatchedInfo<ItemType::node>().m_new_owner, node_owner_vector);
+ return node_owner_vector;
+ }());
this->_buildItemToSubItemDescriptor<NodeOfCell>();
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 73ebc548d98a01e4999611927272206df4af1c23..9ff5913b60799148f99bb6fdf4cff2cedd56973c 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -180,11 +180,11 @@ class ConnectivityDispatcher
template <typename DataType, ItemType item_type, typename ConnectivityPtr>
void _gatherFrom(const ItemValue<DataType, item_type, ConnectivityPtr>& data_to_gather,
- std::vector<std::remove_const_t<DataType>>& gathered_vector);
+ Array<std::remove_const_t<DataType>>& gathered_array);
template <typename DataType, typename ItemOfItem, typename ConnectivityPtr>
void _gatherFrom(const SubItemValuePerItem<DataType, ItemOfItem, ConnectivityPtr>& data_to_gather,
- std::vector<Array<std::remove_const_t<DataType>>>& gathered_vector);
+ Array<std::remove_const_t<DataType>>& gathered_array);
template <typename SubItemOfItemT>
void _buildNumberOfSubItemPerItemToRecvByProc();
diff --git a/src/mesh/ConnectivityMatrix.hpp b/src/mesh/ConnectivityMatrix.hpp
index 6b777288cd6b8e3793a9addff6629c661d45d82d..eadc85347cb756932616e655d7b625a2a15505b7 100644
--- a/src/mesh/ConnectivityMatrix.hpp
+++ b/src/mesh/ConnectivityMatrix.hpp
@@ -70,34 +70,22 @@ class ConnectivityMatrix
#endif // NDEBUG
}
- PUGS_INLINE
- ConnectivityMatrix(const std::vector<std::vector<unsigned int>>& initializer) noexcept : m_is_built{true}
- {
- m_row_map = [&] {
- Array<uint32_t> row_map(initializer.size() + 1);
- row_map[0] = 0;
- for (size_t i = 0; i < initializer.size(); ++i) {
- row_map[i + 1] = row_map[i] + initializer[i].size();
- }
- return row_map;
- }();
-
- m_column_indices = [&] {
- Array<uint32_t> column_indices(m_row_map[m_row_map.size() - 1]);
- size_t index = 0;
- for (const auto& row : initializer) {
- for (const auto& col_index : row) {
- column_indices[index++] = col_index;
- }
- }
- return column_indices;
- }();
- }
-
ConnectivityMatrix& operator=(const ConnectivityMatrix&) = default;
- ConnectivityMatrix& operator=(ConnectivityMatrix&&) = default;
+ ConnectivityMatrix& operator=(ConnectivityMatrix&&) = default;
- ConnectivityMatrix() = default;
+ ConnectivityMatrix(bool is_built = false) : m_is_built{is_built}
+ {
+ // this is useful to build
+ if (is_built) {
+ m_row_map = [&] {
+ Array<uint32_t> row_map(1);
+ row_map[0] = 0;
+ return row_map;
+ }();
+
+ m_column_indices = Array<uint32_t>(0);
+ }
+ }
ConnectivityMatrix(const ConnectivityMatrix&) = default;
ConnectivityMatrix(ConnectivityMatrix&&) = default;
~ConnectivityMatrix() = default;
diff --git a/src/mesh/ConnectivityUtils.cpp b/src/mesh/ConnectivityUtils.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..af0220ac448b76fb2e0101c278719cc07da989bc
--- /dev/null
+++ b/src/mesh/ConnectivityUtils.cpp
@@ -0,0 +1,88 @@
+#include <mesh/ConnectivityUtils.hpp>
+
+#include <mesh/Connectivity.hpp>
+#include <utils/Messenger.hpp>
+
+template <size_t Dimension, ItemType SourceItemType, ItemType TargetItemType>
+bool
+checkItemToHigherDimensionItem(const Connectivity<Dimension>& connectivity)
+{
+ static_assert(SourceItemType < TargetItemType);
+ bool is_valid = true;
+
+ const auto& item_to_item_matrix = connectivity.template getItemToItemMatrix<SourceItemType, TargetItemType>();
+ const auto& item_number = connectivity.template number<TargetItemType>();
+
+ for (ItemIdT<SourceItemType> source_item_id = 0; source_item_id < connectivity.template numberOf<SourceItemType>();
+ ++source_item_id) {
+ auto target_item_list = item_to_item_matrix[source_item_id];
+ for (size_t i = 0; i < target_item_list.size() - 1; ++i) {
+ is_valid &= (item_number[target_item_list[i + 1]] > item_number[target_item_list[i]]);
+ }
+ }
+
+ return is_valid;
+}
+
+template <size_t Dimension>
+bool checkConnectivityOrdering(const Connectivity<Dimension>&);
+
+template <>
+bool
+checkConnectivityOrdering(const Connectivity<1>& connectivity)
+{
+ bool is_valid = true;
+ is_valid &= checkItemToHigherDimensionItem<1, ItemType::node, ItemType::cell>(connectivity);
+
+ return is_valid;
+}
+template <>
+
+bool
+checkConnectivityOrdering(const Connectivity<2>& connectivity)
+{
+ bool is_valid = true;
+ is_valid &= checkItemToHigherDimensionItem<2, ItemType::node, ItemType::face>(connectivity);
+ is_valid &= checkItemToHigherDimensionItem<2, ItemType::node, ItemType::cell>(connectivity);
+
+ is_valid &= checkItemToHigherDimensionItem<2, ItemType::face, ItemType::cell>(connectivity);
+
+ return is_valid;
+}
+
+bool
+checkConnectivityOrdering(const Connectivity<3>& connectivity)
+{
+ bool is_valid = true;
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::edge>(connectivity);
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::face>(connectivity);
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::node, ItemType::cell>(connectivity);
+
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::edge, ItemType::face>(connectivity);
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::edge, ItemType::cell>(connectivity);
+
+ is_valid &= checkItemToHigherDimensionItem<3, ItemType::face, ItemType::cell>(connectivity);
+
+ return parallel::allReduceAnd(is_valid);
+}
+
+bool
+checkConnectivityOrdering(const IConnectivity& connecticity)
+{
+ switch (connecticity.dimension()) {
+ case 1: {
+ return checkConnectivityOrdering(dynamic_cast<const Connectivity<1>&>(connecticity));
+ }
+ case 2: {
+ return checkConnectivityOrdering(dynamic_cast<const Connectivity<2>&>(connecticity));
+ }
+ case 3: {
+ return checkConnectivityOrdering(dynamic_cast<const Connectivity<3>&>(connecticity));
+ }
+ // LCOV_EXCL_START
+ default: {
+ throw UnexpectedError("invalid dimension");
+ }
+ // LCOV_EXCL_STOP
+ }
+}
diff --git a/src/mesh/ConnectivityUtils.hpp b/src/mesh/ConnectivityUtils.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..fa359a7e1c532e3b83c91d3e42d1bd8e3d7a45d8
--- /dev/null
+++ b/src/mesh/ConnectivityUtils.hpp
@@ -0,0 +1,8 @@
+#ifndef CONNECTIVITY_UTILS_HPP
+#define CONNECTIVITY_UTILS_HPP
+
+#include "mesh/IConnectivity.hpp"
+
+bool checkConnectivityOrdering(const IConnectivity&);
+
+#endif // CONNECTIVITY_UTILS_HPP
diff --git a/src/mesh/DiamondDualConnectivityBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
index 71e316461d94a0891b721a0fff85e85c145ef26f..c6a8fd78a9774dc63cbc5e5398ce721c5cfb9acb 100644
--- a/src/mesh/DiamondDualConnectivityBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -11,6 +11,7 @@
#include <utils/Messenger.hpp>
#include <utils/Stringify.hpp>
+#include <optional>
#include <vector>
template <size_t Dimension>
@@ -44,21 +45,21 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
}
}
- diamond_descriptor.node_number_vector.resize(diamond_number_of_nodes);
+ Array<int> node_number_vector(diamond_number_of_nodes);
parallel_for(m_primal_node_to_dual_node_map.size(), [&](size_t i) {
const auto [primal_node_id, diamond_dual_node_id] = m_primal_node_to_dual_node_map[i];
- diamond_descriptor.node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
+ node_number_vector[diamond_dual_node_id] = primal_node_number[primal_node_id];
});
const size_t cell_number_shift = max(primal_node_number) + 1;
parallel_for(primal_number_of_cells, [&](size_t i) {
const auto [primal_cell_id, diamond_dual_node_id] = m_primal_cell_to_dual_node_map[i];
- diamond_descriptor.node_number_vector[diamond_dual_node_id] =
- primal_cell_number[primal_cell_id] + cell_number_shift;
+ node_number_vector[diamond_dual_node_id] = primal_cell_number[primal_cell_id] + cell_number_shift;
});
+ diamond_descriptor.setNodeNumberVector(node_number_vector);
{
m_primal_face_to_dual_cell_map = FaceIdToCellIdMap{primal_number_of_faces};
@@ -68,14 +69,17 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
}
}
- diamond_descriptor.cell_number_vector.resize(diamond_number_of_cells);
- const auto& primal_face_number = primal_connectivity.faceNumber();
- parallel_for(diamond_number_of_cells, [&](size_t i) {
- const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
- diamond_descriptor.cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
- });
+ diamond_descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector(diamond_number_of_cells);
+ const auto& primal_face_number = primal_connectivity.faceNumber();
+ parallel_for(diamond_number_of_cells, [&](size_t i) {
+ const auto [primal_face_id, dual_cell_id] = m_primal_face_to_dual_cell_map[i];
+ cell_number_vector[dual_cell_id] = primal_face_number[primal_face_id];
+ });
+ return cell_number_vector;
+ }());
- diamond_descriptor.cell_type_vector.resize(diamond_number_of_cells);
+ Array<CellType> cell_type_vector(diamond_number_of_cells);
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
@@ -87,97 +91,108 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connec
if constexpr (Dimension == 2) {
if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Triangle;
+ cell_type_vector[i_cell] = CellType::Triangle;
} else {
Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Quadrangle;
+ cell_type_vector[i_cell] = CellType::Quadrangle;
}
} else if constexpr (Dimension == 3) {
if (primal_face_cell_list.size() == 1) {
if (primal_face_to_node_matrix[face_id].size() == 3) {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Tetrahedron;
+ cell_type_vector[i_cell] = CellType::Tetrahedron;
} else {
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Pyramid;
+ cell_type_vector[i_cell] = CellType::Pyramid;
}
} else {
Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_type_vector[i_cell] = CellType::Diamond;
+ cell_type_vector[i_cell] = CellType::Diamond;
}
}
});
- diamond_descriptor.cell_to_node_vector.resize(diamond_number_of_cells);
+ diamond_descriptor.setCellTypeVector(cell_type_vector);
- const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
- const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
- parallel_for(primal_number_of_faces, [&](FaceId face_id) {
- const size_t& i_diamond_cell = face_id;
- const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
- const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
- if (primal_face_cell_list.size() == 1) {
- diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 1);
-
- const CellId cell_id = primal_face_cell_list[0];
- const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
+ Array<const unsigned int> cell_to_node_row = [&] {
+ Array<unsigned int> tmp_cell_to_node_row(primal_number_of_faces + 1);
+ tmp_cell_to_node_row[0] = 0;
+ for (FaceId face_id = 0; face_id < primal_number_of_faces; ++face_id) {
+ tmp_cell_to_node_row[face_id + 1] = tmp_cell_to_node_row[face_id] + primal_face_to_node_matrix[face_id].size() +
+ primal_face_to_cell_matrix[face_id].size();
+ }
- for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node] = primal_face_node_list[i_node];
- }
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size()] =
- primal_number_of_nodes + cell_id;
+ return tmp_cell_to_node_row;
+ }();
+
+ Array<const unsigned int> cell_to_node_list = [&] {
+ Array<unsigned int> tmp_cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+ const auto& primal_face_local_number_in_their_cells = primal_connectivity.faceLocalNumbersInTheirCells();
+ const auto& cell_face_is_reversed = primal_connectivity.cellFaceIsReversed();
+ parallel_for(primal_number_of_faces, [&](FaceId face_id) {
+ const auto& primal_face_cell_list = primal_face_to_cell_matrix[face_id];
+ const auto& primal_face_node_list = primal_face_to_node_matrix[face_id];
+ const size_t first_node = cell_to_node_row[face_id];
+ if (primal_face_cell_list.size() == 1) {
+ const CellId cell_id = primal_face_cell_list[0];
+ const auto i_face_in_cell = primal_face_local_number_in_their_cells(face_id, 0);
- if constexpr (Dimension == 2) {
- if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][1]);
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ tmp_cell_to_node_list[first_node + i_node] = primal_face_node_list[i_node];
}
- } else {
- if (not cell_face_is_reversed(cell_id, i_face_in_cell)) {
- // In 3D the basis of the pyramid is described in the
- // indirect way IF the face is not reversed. In other words
- // the "topological normal" must point to the "top" of the
- // pyramid.
- for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node],
- diamond_descriptor
- .cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() - 1 - i_node]);
+ tmp_cell_to_node_list[first_node + primal_face_node_list.size()] = primal_number_of_nodes + cell_id;
+
+ if constexpr (Dimension == 2) {
+ if (cell_face_is_reversed(cell_id, i_face_in_cell)) {
+ std::swap(tmp_cell_to_node_list[first_node], tmp_cell_to_node_list[first_node + 1]);
+ }
+ } else {
+ if (not cell_face_is_reversed(cell_id, i_face_in_cell)) {
+ // In 3D the basis of the pyramid is described in the
+ // indirect way IF the face is not reversed. In other words
+ // the "topological normal" must point to the "top" of the
+ // pyramid.
+ for (size_t i_node = 0; i_node < primal_face_node_list.size() / 2; ++i_node) {
+ std::swap(tmp_cell_to_node_list[first_node + i_node],
+ tmp_cell_to_node_list[first_node + primal_face_node_list.size() - 1 - i_node]);
+ }
}
- }
- }
- } else {
- Assert(primal_face_cell_list.size() == 2);
- diamond_descriptor.cell_to_node_vector[i_diamond_cell].resize(primal_face_node_list.size() + 2);
-
- const CellId cell0_id = primal_face_cell_list[0];
- const CellId cell1_id = primal_face_cell_list[1];
- const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
-
- if constexpr (Dimension == 2) {
- Assert(primal_face_node_list.size() == 2);
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][1] = primal_face_node_list[0];
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][2] = primal_number_of_nodes + cell1_id;
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][3] = primal_face_node_list[1];
-
- if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][1],
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][3]);
}
} else {
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][0] = primal_number_of_nodes + cell0_id;
- for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][i_node + 1] = primal_face_node_list[i_node];
- }
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1] =
- primal_number_of_nodes + cell1_id;
+ Assert(primal_face_cell_list.size() == 2);
- if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
- std::swap(diamond_descriptor.cell_to_node_vector[i_diamond_cell][0],
- diamond_descriptor.cell_to_node_vector[i_diamond_cell][primal_face_node_list.size() + 1]);
+ const CellId cell0_id = primal_face_cell_list[0];
+ const CellId cell1_id = primal_face_cell_list[1];
+ const auto i_face_in_cell0 = primal_face_local_number_in_their_cells(face_id, 0);
+
+ if constexpr (Dimension == 2) {
+ Assert(primal_face_node_list.size() == 2);
+
+ tmp_cell_to_node_list[first_node + 0] = primal_number_of_nodes + cell0_id;
+ tmp_cell_to_node_list[first_node + 1] = primal_face_node_list[0];
+ tmp_cell_to_node_list[first_node + 2] = primal_number_of_nodes + cell1_id;
+ tmp_cell_to_node_list[first_node + 3] = primal_face_node_list[1];
+
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+ std::swap(tmp_cell_to_node_list[first_node + 1], tmp_cell_to_node_list[first_node + 3]);
+ }
+ } else {
+ tmp_cell_to_node_list[first_node + 0] = primal_number_of_nodes + cell0_id;
+ for (size_t i_node = 0; i_node < primal_face_node_list.size(); ++i_node) {
+ tmp_cell_to_node_list[first_node + i_node + 1] = primal_face_node_list[i_node];
+ }
+ tmp_cell_to_node_list[first_node + primal_face_node_list.size() + 1] = primal_number_of_nodes + cell1_id;
+
+ if (cell_face_is_reversed(cell0_id, i_face_in_cell0)) {
+ std::swap(tmp_cell_to_node_list[first_node],
+ tmp_cell_to_node_list[first_node + primal_face_node_list.size() + 1]);
+ }
}
}
- }
- });
+ });
+
+ return tmp_cell_to_node_list;
+ }();
+
+ diamond_descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
}
template <size_t Dimension>
@@ -198,11 +213,12 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
}
+ const auto& node_number_vector = diamond_descriptor.nodeNumberVector();
{
const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
std::unordered_map<unsigned int, NodeId> node_to_id_map;
- for (size_t i_node = 0; i_node < diamond_descriptor.node_number_vector.size(); ++i_node) {
- node_to_id_map[diamond_descriptor.node_number_vector[i_node]] = i_node;
+ for (size_t i_node = 0; i_node < node_number_vector.size(); ++i_node) {
+ node_to_id_map[node_number_vector[i_node]] = i_node;
}
return node_to_id_map;
}();
@@ -238,20 +254,32 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
}
{
- const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
-
- using Face = ConnectivityFace<Dimension>;
-
- 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 < diamond_descriptor.face_to_node_vector.size(); ++l) {
- const auto& node_vector = diamond_descriptor.face_to_node_vector[l];
-
- face_to_id_map[Face(node_vector, diamond_descriptor.node_number_vector)] = l;
+ const auto& primal_face_to_node_matrix = primal_connectivity.faceToNodeMatrix();
+ const auto& diamond_node_to_face_matrix = diamond_descriptor.nodeToFaceMatrix();
+ const auto& diamond_face_to_node_matrix = diamond_descriptor.faceToNodeMatrix();
+
+ const auto find_face = [&](std::vector<uint32_t> node_list) -> std::optional<FaceId> {
+ // The node list of already sorted correctly
+ const auto& face_id_vector = diamond_node_to_face_matrix[node_list[0]];
+
+ for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+ const FaceId face_id = face_id_vector[i_face];
+ const auto& face_node_id_list = diamond_face_to_node_matrix[face_id];
+ if (face_node_id_list.size() == node_list.size()) {
+ bool is_same = true;
+ for (size_t i_node = 1; i_node < face_node_id_list.size(); ++i_node) {
+ is_same &= (face_node_id_list[i_node] == node_list[i_node]);
+ }
+ if (is_same) {
+ return face_id;
+ }
+ }
}
- return face_to_id_map;
- }();
+ return std::nullopt;
+ };
+
+ std::vector<unsigned int> face_node_list;
for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
++i_face_list) {
const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
@@ -265,16 +293,15 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
const auto& primal_face_node_list = primal_face_to_node_matrix[primal_face_id];
- const auto i_diamond_face = [&]() {
- std::vector<unsigned int> node_list(primal_face_node_list.size());
- for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
- node_list[i] = primal_face_node_list[i];
- }
- return face_to_id_map.find(Face(node_list, diamond_descriptor.node_number_vector));
- }();
+ face_node_list.clear();
+
+ for (size_t i = 0; i < primal_face_node_list.size(); ++i) {
+ face_node_list.push_back(primal_face_node_list[i]);
+ }
- if (i_diamond_face != face_to_id_map.end()) {
- diamond_face_list.push_back(i_diamond_face->second);
+ auto face_id = find_face(face_node_list);
+ if (face_id.has_value()) {
+ diamond_face_list.push_back(face_id.value());
}
}
return diamond_face_list;
@@ -293,29 +320,40 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
if constexpr (Dimension == 3) {
const auto& primal_edge_to_node_matrix = primal_connectivity.edgeToNodeMatrix();
- using Edge = ConnectivityFace<2>;
-
- 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 < diamond_descriptor.edge_to_node_vector.size(); ++l) {
- const auto& node_vector = diamond_descriptor.edge_to_node_vector[l];
- edge_to_id_map[Edge(node_vector, diamond_descriptor.node_number_vector)] = l;
- }
- return edge_to_id_map;
- }();
{
- const size_t number_of_edges = diamond_descriptor.edge_to_node_vector.size();
- diamond_descriptor.edge_number_vector.resize(number_of_edges);
- for (size_t i_edge = 0; i_edge < number_of_edges; ++i_edge) {
- diamond_descriptor.edge_number_vector[i_edge] = i_edge;
- }
+ const size_t number_of_edges = diamond_descriptor.edgeToNodeMatrix().numberOfRows();
+ diamond_descriptor.setEdgeNumberVector([&] {
+ Array<int> edge_number_vector(number_of_edges);
+ parallel_for(
+ number_of_edges, PUGS_LAMBDA(size_t i_edge) { edge_number_vector[i_edge] = i_edge; });
+ return edge_number_vector;
+ }());
+
// LCOV_EXCL_START
if (parallel::size() > 1) {
throw NotImplementedError("parallel edge numbering is undefined");
}
// LCOV_EXCL_STOP
}
+ const auto& diamond_node_to_edge_matrix = diamond_descriptor.nodeToEdgeMatrix();
+ const auto& diamond_edge_to_node_matrix = diamond_descriptor.edgeToNodeMatrix();
+
+ const auto find_edge = [&](uint32_t node0, uint32_t node1) -> std::optional<EdgeId> {
+ if (node_number_vector[node0] > node_number_vector[node1]) {
+ std::swap(node0, node1);
+ }
+ const auto& edge_id_vector = diamond_node_to_edge_matrix[node0];
+
+ for (size_t i_edge = 0; i_edge < edge_id_vector.size(); ++i_edge) {
+ const EdgeId edge_id = edge_id_vector[i_edge];
+ if (diamond_edge_to_node_matrix[edge_id][1] == node1) {
+ return edge_id;
+ }
+ }
+
+ return std::nullopt;
+ };
for (size_t i_edge_list = 0; i_edge_list < primal_connectivity.template numberOfRefItemList<ItemType::edge>();
++i_edge_list) {
@@ -330,16 +368,10 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
const auto& primal_edge_node_list = primal_edge_to_node_matrix[primal_edge_id];
- const auto i_diamond_edge = [&]() {
- std::vector<unsigned int> node_list(primal_edge_node_list.size());
- for (size_t i = 0; i < primal_edge_node_list.size(); ++i) {
- node_list[i] = primal_edge_node_list[i];
- }
- return edge_to_id_map.find(Edge(node_list, diamond_descriptor.node_number_vector));
- }();
+ const auto diamond_edge_id = find_edge(primal_edge_node_list[0], primal_edge_node_list[1]);
- if (i_diamond_edge != edge_to_id_map.end()) {
- diamond_edge_list.push_back(i_diamond_edge->second);
+ if (diamond_edge_id.has_value()) {
+ diamond_edge_list.push_back(diamond_edge_id.value());
}
}
return diamond_edge_list;
@@ -359,62 +391,64 @@ DiamondDualConnectivityBuilder::_buildDiamondConnectivityFrom(const IConnectivit
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
- diamond_descriptor.node_owner_vector.resize(diamond_descriptor.node_number_vector.size());
-
- {
+ diamond_descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(node_number_vector.size());
const auto& primal_node_owner = primal_connectivity.nodeOwner();
for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
- diamond_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
+ node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
}
const auto& primal_cell_owner = primal_connectivity.cellOwner();
for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- diamond_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+ node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
}
- }
+ return node_owner_vector;
+ }());
- {
- diamond_descriptor.cell_owner_vector.resize(diamond_descriptor.cell_number_vector.size());
+ diamond_descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(diamond_descriptor.cellNumberVector().size());
const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
const auto& primal_face_owner = primal_connectivity.faceOwner();
- for (FaceId primal_face_id = 0; primal_face_id < primal_number_of_faces; ++primal_face_id) {
- diamond_descriptor.cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
- }
- }
-
- {
- std::vector<int> face_cell_owner(diamond_descriptor.face_number_vector.size());
- std::fill(std::begin(face_cell_owner), std::end(face_cell_owner), parallel::size());
-
- for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
- const auto& cell_face_list = diamond_descriptor.cell_to_face_vector[i_cell];
+ parallel_for(
+ primal_number_of_faces, PUGS_LAMBDA(const FaceId primal_face_id) {
+ cell_owner_vector[primal_face_id] = primal_face_owner[primal_face_id];
+ });
+ return cell_owner_vector;
+ }());
+
+ const auto& diamond_cell_owner_vector = diamond_descriptor.cellOwnerVector();
+ const auto& diamond_cell_to_face_matrix = diamond_descriptor.cellToFaceMatrix();
+
+ diamond_descriptor.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector(diamond_descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::rank());
+
+ for (size_t i_cell = 0; i_cell < diamond_cell_to_face_matrix.numberOfRows(); ++i_cell) {
+ const auto& cell_face_list = diamond_cell_to_face_matrix[i_cell];
for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- const size_t face_id = cell_face_list[i_face];
- face_cell_owner[face_id] = std::min(face_cell_owner[face_id], diamond_descriptor.cell_number_vector[i_cell]);
+ const size_t face_id = cell_face_list[i_face];
+ face_owner_vector[face_id] = std::min(face_owner_vector[face_id], diamond_cell_owner_vector[i_cell]);
}
}
-
- diamond_descriptor.face_owner_vector.resize(face_cell_owner.size());
- for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) {
- diamond_descriptor.face_owner_vector[i_face] = diamond_descriptor.cell_owner_vector[face_cell_owner[i_face]];
- }
- }
+ return face_owner_vector;
+ }());
if constexpr (Dimension == 3) {
- std::vector<int> edge_cell_owner(diamond_descriptor.edge_number_vector.size());
- std::fill(std::begin(edge_cell_owner), std::end(edge_cell_owner), parallel::size());
-
- for (size_t i_cell = 0; i_cell < diamond_descriptor.cell_to_face_vector.size(); ++i_cell) {
- const auto& cell_edge_list = diamond_descriptor.cell_to_edge_vector[i_cell];
- for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
- const size_t edge_id = cell_edge_list[i_edge];
- edge_cell_owner[edge_id] = std::min(edge_cell_owner[edge_id], diamond_descriptor.cell_number_vector[i_cell]);
+ const auto& diamond_cell_to_edge_matrix = diamond_descriptor.cellToEdgeMatrix();
+
+ diamond_descriptor.setEdgeOwnerVector([&] {
+ Array<int> edge_owner_vector(diamond_descriptor.edgeNumberVector().size());
+ edge_owner_vector.fill(parallel::rank());
+
+ for (size_t i_cell = 0; i_cell < diamond_cell_to_edge_matrix.numberOfRows(); ++i_cell) {
+ const auto& cell_edge_list = diamond_cell_to_edge_matrix[i_cell];
+ for (size_t i_edge = 0; i_edge < cell_edge_list.size(); ++i_edge) {
+ const size_t edge_id = cell_edge_list[i_edge];
+ edge_owner_vector[edge_id] = std::min(edge_owner_vector[edge_id], diamond_cell_owner_vector[i_cell]);
+ }
}
- }
- diamond_descriptor.edge_owner_vector.resize(edge_cell_owner.size());
- for (size_t i_edge = 0; i_edge < edge_cell_owner.size(); ++i_edge) {
- diamond_descriptor.edge_owner_vector[i_edge] = diamond_descriptor.cell_owner_vector[edge_cell_owner[i_edge]];
- }
+ return edge_owner_vector;
+ }());
}
m_connectivity = ConnectivityType::build(diamond_descriptor);
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualMeshBuilder.cpp
index 21e413ff4e94519b79fbc26e87d742e67c0436c2..64711edc72055effd9e5608bc01b71ea184ee1fc 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualMeshBuilder.cpp
@@ -43,8 +43,6 @@ DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const IMesh& i_mesh)
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const IMesh& i_mesh)
{
- std::cout << "building DiamondDualMesh\n";
-
switch (i_mesh.dimension()) {
case 2: {
this->_buildDualDiamondMeshFrom<2>(i_mesh);
diff --git a/src/mesh/Dual1DConnectivityBuilder.cpp b/src/mesh/Dual1DConnectivityBuilder.cpp
index 30c465b8b07375fce9f82fb76a26ac24848f3c67..52f3ddc977adc531496bf02b10e2b3ef36dfa5f1 100644
--- a/src/mesh/Dual1DConnectivityBuilder.cpp
+++ b/src/mesh/Dual1DConnectivityBuilder.cpp
@@ -28,14 +28,14 @@ Dual1DConnectivityBuilder::_buildConnectivityDescriptor(const Connectivity<1>& p
const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
size_t next_kept_node_id = 0;
- dual_descriptor.node_number_vector.resize(dual_number_of_nodes);
+ Array<int> node_number_vector(dual_number_of_nodes);
const auto& primal_node_number = primal_connectivity.nodeNumber();
for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
const auto& primal_node_cell_list = primal_node_to_cell_matrix[primal_node_id];
if (primal_node_cell_list.size() == 1) {
- dual_descriptor.node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id];
+ node_number_vector[next_kept_node_id++] = primal_node_number[primal_node_id];
}
}
@@ -46,43 +46,56 @@ Dual1DConnectivityBuilder::_buildConnectivityDescriptor(const Connectivity<1>& p
const size_t cell_number_shift = max(primal_node_number) + 1;
for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- dual_descriptor.node_number_vector[primal_number_of_kept_nodes + primal_cell_id] =
+ node_number_vector[primal_number_of_kept_nodes + primal_cell_id] =
primal_cell_number[primal_cell_id] + cell_number_shift;
}
+ dual_descriptor.setNodeNumberVector(node_number_vector);
Assert(number_of_kept_nodes == next_kept_node_id, "unexpected number of kept nodes");
- dual_descriptor.cell_number_vector.resize(dual_number_of_cells);
- for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
- dual_descriptor.cell_number_vector[primal_node_id] = primal_node_number[primal_node_id];
- }
-
- dual_descriptor.cell_type_vector = std::vector<CellType>(dual_number_of_cells, CellType::Line);
+ dual_descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector(dual_number_of_cells);
+ parallel_for(
+ primal_number_of_nodes, PUGS_LAMBDA(const NodeId primal_node_id) {
+ cell_number_vector[primal_node_id] = primal_node_number[primal_node_id];
+ });
+ return cell_number_vector;
+ }());
- dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
+ Array<CellType> cell_type_vector(dual_number_of_cells);
+ cell_type_vector.fill(CellType::Line);
+ dual_descriptor.setCellTypeVector(cell_type_vector);
const auto& primal_node_local_number_in_their_cells = primal_connectivity.nodeLocalNumbersInTheirCells();
+ Array<unsigned int> cell_to_node_row(dual_number_of_cells + 1);
+ parallel_for(
+ cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
+
+ Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
{
size_t next_kept_node_id = 0;
+ size_t i_cell_node = 0;
for (NodeId i_node = 0; i_node < primal_connectivity.numberOfNodes(); ++i_node) {
- const size_t& i_dual_cell = i_node;
const auto& primal_node_cell_list = primal_node_to_cell_matrix[i_node];
- dual_descriptor.cell_to_node_vector[i_dual_cell].resize(2);
if (primal_node_cell_list.size() == 1) {
const auto i_node_in_cell = primal_node_local_number_in_their_cells(i_node, 0);
- dual_descriptor.cell_to_node_vector[i_dual_cell][i_node_in_cell] = next_kept_node_id++;
- dual_descriptor.cell_to_node_vector[i_dual_cell][1 - i_node_in_cell] =
- number_of_kept_nodes + primal_node_cell_list[0];
+ cell_to_node_list[i_cell_node + i_node_in_cell] = next_kept_node_id++;
+ cell_to_node_list[i_cell_node + 1 - i_node_in_cell] = number_of_kept_nodes + primal_node_cell_list[0];
+
+ i_cell_node += 2;
+
} else {
const auto i0 = primal_node_local_number_in_their_cells(i_node, 0);
- dual_descriptor.cell_to_node_vector[i_dual_cell][0] = number_of_kept_nodes + primal_node_cell_list[1 - i0];
- dual_descriptor.cell_to_node_vector[i_dual_cell][1] = number_of_kept_nodes + primal_node_cell_list[i0];
+ cell_to_node_list[i_cell_node++] = number_of_kept_nodes + primal_node_cell_list[1 - i0];
+ cell_to_node_list[i_cell_node++] = number_of_kept_nodes + primal_node_cell_list[i0];
}
}
}
+
+ dual_descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
}
void
@@ -97,10 +110,12 @@ Dual1DConnectivityBuilder::_buildConnectivityFrom(const IConnectivity& i_primal_
this->_buildConnectivityDescriptor(primal_connectivity, dual_descriptor);
{
+ const auto& dual_node_number_vector = dual_descriptor.nodeNumberVector();
+
const std::unordered_map<unsigned int, NodeId> node_to_id_map = [&] {
std::unordered_map<unsigned int, NodeId> node_to_id_map;
- for (size_t i_node = 0; i_node < dual_descriptor.node_number_vector.size(); ++i_node) {
- node_to_id_map[dual_descriptor.node_number_vector[i_node]] = i_node;
+ for (size_t i_node = 0; i_node < dual_node_number_vector.size(); ++i_node) {
+ node_to_id_map[dual_node_number_vector[i_node]] = i_node;
}
return node_to_id_map;
}();
@@ -138,31 +153,33 @@ Dual1DConnectivityBuilder::_buildConnectivityFrom(const IConnectivity& i_primal_
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
- dual_descriptor.node_owner_vector.resize(dual_descriptor.node_number_vector.size());
-
- {
+ dual_descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(dual_descriptor.nodeNumberVector().size());
const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
const auto& primal_node_owner = primal_connectivity.nodeOwner();
size_t next_kept_node_id = 0;
for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
if (node_to_cell_matrix[primal_node_id].size() == 1) {
- dual_descriptor.node_owner_vector[next_kept_node_id++] = primal_node_owner[primal_node_id];
+ node_owner_vector[next_kept_node_id++] = primal_node_owner[primal_node_id];
}
}
const size_t number_of_kept_nodes = next_kept_node_id;
const auto& primal_cell_owner = primal_connectivity.cellOwner();
for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- dual_descriptor.node_owner_vector[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
+ node_owner_vector[number_of_kept_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
}
- }
- {
- dual_descriptor.cell_owner_vector.resize(dual_descriptor.cell_number_vector.size());
+ return node_owner_vector;
+ }());
+
+ dual_descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(dual_descriptor.cellNumberVector().size());
const auto& primal_node_owner = primal_connectivity.nodeOwner();
- for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
- dual_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
- }
- }
+ parallel_for(
+ primal_number_of_nodes,
+ PUGS_LAMBDA(NodeId primal_node_id) { cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id]; });
+ return cell_owner_vector;
+ }());
m_connectivity = ConnectivityType::build(dual_descriptor);
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index e555ca224c1d02f11e11696542774f2e470b6432..e4ae59ea7cbb6fa2ec3f622c0868ff79ff9f8bbf 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -35,19 +35,29 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
{
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector = gmsh_data.__verticesNumbers;
- descriptor.cell_type_vector.resize(nb_cells);
- descriptor.cell_number_vector.resize(nb_cells);
- descriptor.cell_to_node_vector.resize(nb_cells);
+ descriptor.setNodeNumberVector(convert_to_array(gmsh_data.__verticesNumbers));
+ Array<CellType> cell_type_vector(nb_cells);
+ Array<int> cell_number_vector(nb_cells);
- for (size_t j = 0; j < nb_cells; ++j) {
- descriptor.cell_to_node_vector[j].resize(2);
- for (int r = 0; r < 2; ++r) {
- descriptor.cell_to_node_vector[j][r] = gmsh_data.__edges[j][r];
+ Array<unsigned int> cell_to_node_row(nb_cells + 1);
+ parallel_for(
+ cell_to_node_row.size(), PUGS_LAMBDA(const CellId cell_id) { cell_to_node_row[cell_id] = 2 * cell_id; });
+
+ Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+ for (CellId cell_id = 0; cell_id < nb_cells; ++cell_id) {
+ for (int i_node = 0; i_node < 2; ++i_node) {
+ cell_to_node_list[2 * cell_id + i_node] = gmsh_data.__edges[cell_id][i_node];
}
- descriptor.cell_type_vector[j] = CellType::Line;
- descriptor.cell_number_vector[j] = gmsh_data.__edges_number[j];
}
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
+
+ cell_type_vector.fill(CellType::Line);
+ descriptor.setCellTypeVector(cell_type_vector);
+
+ for (size_t j = 0; j < nb_cells; ++j) {
+ cell_number_vector[j] = gmsh_data.__edges_number[j];
+ }
+ descriptor.setCellNumberVector(cell_number_vector);
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
@@ -56,12 +66,14 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
ref_points_map[ref].push_back(point_number);
}
- Array<size_t> node_nb_cell(descriptor.node_number_vector.size());
+ Array<size_t> node_nb_cell(descriptor.nodeNumberVector().size());
node_nb_cell.fill(0);
+ const auto& cell_to_node_matrix = descriptor.cellToNodeMatrix();
for (size_t j = 0; j < nb_cells; ++j) {
+ const auto& cell_node_list = cell_to_node_matrix[j];
for (int r = 0; r < 2; ++r) {
- node_nb_cell[descriptor.cell_to_node_vector[j][r]] += 1;
+ node_nb_cell[cell_node_list[r]] += 1;
}
}
@@ -116,12 +128,18 @@ GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData&
descriptor.addRefItemList(RefCellList(ref_id, cell_list, false));
}
- 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());
-
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(nb_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
+
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(descriptor.nodeNumberVector().size());
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
+ ;
m_connectivity = Connectivity1D::build(descriptor);
}
@@ -130,32 +148,56 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
{
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector = gmsh_data.__verticesNumbers;
- descriptor.cell_type_vector.resize(nb_cells);
- descriptor.cell_number_vector.resize(nb_cells);
- descriptor.cell_to_node_vector.resize(nb_cells);
+ descriptor.setNodeNumberVector(convert_to_array(gmsh_data.__verticesNumbers));
+ Array<CellType> cell_type_vector(nb_cells);
+ Array<int> cell_number_vector(nb_cells);
- const size_t nb_triangles = gmsh_data.__triangles.size();
- for (size_t j = 0; j < nb_triangles; ++j) {
- descriptor.cell_to_node_vector[j].resize(3);
- for (int r = 0; r < 3; ++r) {
- descriptor.cell_to_node_vector[j][r] = gmsh_data.__triangles[j][r];
+ Array<unsigned int> cell_to_node_row(gmsh_data.__triangles.size() + gmsh_data.__quadrangles.size() + 1);
+ {
+ cell_to_node_row[0] = 0;
+ size_t i_cell = 0;
+ for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 3;
+ }
+ for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 4;
}
- descriptor.cell_type_vector[j] = CellType::Triangle;
- descriptor.cell_number_vector[j] = gmsh_data.__triangles_number[j];
+ }
+
+ Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+ {
+ size_t i_cell_node = 0;
+ for (size_t i_triangle = 0; i_triangle < gmsh_data.__triangles.size(); ++i_triangle) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__triangles[i_triangle][2];
+ }
+ for (size_t i_quadrangle = 0; i_quadrangle < gmsh_data.__quadrangles.size(); ++i_quadrangle) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][2];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__quadrangles[i_quadrangle][3];
+ }
+ }
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
+
+ const size_t nb_triangles = gmsh_data.__triangles.size();
+ for (size_t i_triangle = 0; i_triangle < nb_triangles; ++i_triangle) {
+ cell_type_vector[i_triangle] = CellType::Triangle;
+ cell_number_vector[i_triangle] = gmsh_data.__triangles_number[i_triangle];
}
const size_t nb_quadrangles = gmsh_data.__quadrangles.size();
- for (size_t j = 0; j < nb_quadrangles; ++j) {
- const size_t jq = j + nb_triangles;
- descriptor.cell_to_node_vector[jq].resize(4);
- for (int r = 0; r < 4; ++r) {
- descriptor.cell_to_node_vector[jq][r] = gmsh_data.__quadrangles[j][r];
- }
- descriptor.cell_type_vector[jq] = CellType::Quadrangle;
- descriptor.cell_number_vector[jq] = gmsh_data.__quadrangles_number[j];
+ for (size_t i_quadrangle = 0; i_quadrangle < nb_quadrangles; ++i_quadrangle) {
+ const size_t i_cell = i_quadrangle + nb_triangles;
+
+ cell_type_vector[i_cell] = CellType::Quadrangle;
+ cell_number_vector[i_cell] = gmsh_data.__quadrangles_number[i_quadrangle];
}
+ descriptor.setCellNumberVector(cell_number_vector);
+ descriptor.setCellTypeVector(cell_type_vector);
+
std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
for (unsigned int j = 0; j < gmsh_data.__triangles_ref.size(); ++j) {
const unsigned int elem_number = j;
@@ -189,64 +231,73 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(descriptor);
- 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;
- 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 = [&] {
+ const auto& face_number_vector = descriptor.faceNumberVector();
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;
+ for (size_t l = 0; l < face_number_vector.size(); ++l) {
+ face_number_id_map[face_number_vector[l]] = l;
}
- Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
+ Assert(face_number_id_map.size() == face_number_vector.size());
return face_number_id_map;
}();
+ const auto& node_number_vector = descriptor.nodeNumberVector();
+ const auto& node_to_face_matrix = descriptor.nodeToFaceMatrix();
+ const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
+
+ const auto find_face = [&](uint32_t node0, uint32_t node1) {
+ if (node_number_vector[node0] > node_number_vector[node1]) {
+ std::swap(node0, node1);
+ }
+ const auto& face_id_vector = node_to_face_matrix[node0];
+
+ for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+ const FaceId face_id = face_id_vector[i_face];
+ if (face_to_node_matrix[face_id][1] == node1) {
+ return face_id;
+ }
+ }
+
+ std::stringstream error_msg;
+ error_msg << "face (" << node0 << ',' << node1 << ") not found";
+ throw NormalError(error_msg.str());
+ };
+
+ Array<int> face_number_vector = copy(descriptor.faceNumberVector());
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
- const unsigned int edge_id = [&] {
- auto i = face_to_id_map.find(Face({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
- if (i == face_to_id_map.end()) {
- std::stringstream error_msg;
- error_msg << "face " << gmsh_data.__edges[e][0] << " not found";
- throw NormalError(error_msg.str());
- }
- return i->second;
- }();
+ const unsigned int edge_id = find_face(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+
const unsigned int& ref = gmsh_data.__edges_ref[e];
ref_faces_map[ref].push_back(edge_id);
- if (descriptor.face_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
+ if (face_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
if (auto i_face = face_number_id_map.find(gmsh_data.__edges_number[e]); i_face != face_number_id_map.end()) {
const int other_edge_id = i_face->second;
- std::swap(descriptor.face_number_vector[edge_id], descriptor.face_number_vector[other_edge_id]);
+ std::swap(face_number_vector[edge_id], face_number_vector[other_edge_id]);
- face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
- face_number_id_map.erase(descriptor.face_number_vector[other_edge_id]);
+ face_number_id_map.erase(face_number_vector[edge_id]);
+ face_number_id_map.erase(face_number_vector[other_edge_id]);
- face_number_id_map[descriptor.face_number_vector[edge_id]] = edge_id;
- face_number_id_map[descriptor.face_number_vector[other_edge_id]] = other_edge_id;
+ face_number_id_map[face_number_vector[edge_id]] = edge_id;
+ face_number_id_map[face_number_vector[other_edge_id]] = other_edge_id;
} else {
- face_number_id_map.erase(descriptor.face_number_vector[edge_id]);
- descriptor.face_number_vector[edge_id] = gmsh_data.__edges_number[e];
- face_number_id_map[descriptor.face_number_vector[edge_id]] = edge_id;
+ face_number_id_map.erase(face_number_vector[edge_id]);
+ face_number_vector[edge_id] = gmsh_data.__edges_number[e];
+ face_number_id_map[face_number_vector[edge_id]] = edge_id;
}
}
}
+ descriptor.setFaceNumberVector(face_number_vector);
- Array<size_t> face_nb_cell(descriptor.face_number_vector.size());
+ Array<size_t> face_nb_cell(descriptor.faceNumberVector().size());
face_nb_cell.fill(0);
- for (size_t j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
- for (size_t l = 0; l < descriptor.cell_to_face_vector[j].size(); ++l) {
- face_nb_cell[descriptor.cell_to_face_vector[j][l]] += 1;
+ const auto& cell_to_face_matrix = descriptor.cellToFaceMatrix();
+ for (size_t j = 0; j < cell_to_face_matrix.numberOfRows(); ++j) {
+ const auto& cell_face_list = cell_to_face_matrix[j];
+ for (size_t l = 0; l < cell_face_list.size(); ++l) {
+ face_nb_cell[cell_face_list[l]] += 1;
}
}
@@ -276,13 +327,14 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
descriptor.addRefItemList(RefFaceList{ref_id, face_list, is_boundary});
}
- Array<bool> is_boundary_node(descriptor.node_number_vector.size());
+ Array<bool> is_boundary_node(node_number_vector.size());
is_boundary_node.fill(false);
for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
if (face_nb_cell[i_face] == 1) {
- for (size_t node_id : descriptor.face_to_node_vector[i_face]) {
- is_boundary_node[node_id] = true;
+ const auto& face_node_list = face_to_node_matrix[i_face];
+ for (size_t i_node = 0; i_node < face_node_list.size(); ++i_node) {
+ is_boundary_node[face_node_list[i_node]] = true;
}
}
}
@@ -319,14 +371,23 @@ GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData&
descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
}
- descriptor.cell_owner_vector.resize(nb_cells);
- std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(nb_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
- 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.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector(descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::rank());
+ return face_owner_vector;
+ }());
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(node_number_vector.size());
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
m_connectivity = Connectivity2D::build(descriptor);
}
@@ -336,54 +397,99 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
{
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector = gmsh_data.__verticesNumbers;
- descriptor.cell_type_vector.resize(nb_cells);
- descriptor.cell_number_vector.resize(nb_cells);
- descriptor.cell_to_node_vector.resize(nb_cells);
+ descriptor.setNodeNumberVector(convert_to_array(gmsh_data.__verticesNumbers));
+ Array<CellType> cell_type_vector(nb_cells);
+ Array<int> cell_number_vector(nb_cells);
const size_t nb_tetrahedra = gmsh_data.__tetrahedra.size();
- for (size_t j = 0; j < nb_tetrahedra; ++j) {
- descriptor.cell_to_node_vector[j].resize(4);
- for (int r = 0; r < 4; ++r) {
- descriptor.cell_to_node_vector[j][r] = gmsh_data.__tetrahedra[j][r];
+ const size_t nb_hexahedra = gmsh_data.__hexahedra.size();
+ const size_t nb_prisms = gmsh_data.__prisms.size();
+ const size_t nb_pyramids = gmsh_data.__pyramids.size();
+
+ Array<unsigned int> cell_to_node_row(nb_cells + 1);
+ {
+ cell_to_node_row[0] = 0;
+ size_t i_cell = 0;
+ for (size_t i_tetrahedron = 0; i_tetrahedron < nb_tetrahedra; ++i_tetrahedron, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 4;
+ }
+ for (size_t i_hexahedron = 0; i_hexahedron < nb_hexahedra; ++i_hexahedron, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 8;
+ }
+ for (size_t i_prism = 0; i_prism < nb_prisms; ++i_prism, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 6;
+ }
+ for (size_t i_pyramid = 0; i_pyramid < nb_pyramids; ++i_pyramid, ++i_cell) {
+ cell_to_node_row[i_cell + 1] = cell_to_node_row[i_cell] + 5;
+ }
+ }
+
+ Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
+ {
+ size_t i_cell_node = 0;
+ for (size_t i_tetrahedron = 0; i_tetrahedron < nb_tetrahedra; ++i_tetrahedron) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][2];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__tetrahedra[i_tetrahedron][3];
+ }
+ for (size_t i_hexahedron = 0; i_hexahedron < nb_hexahedra; ++i_hexahedron) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][2];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][3];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][4];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][5];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][6];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__hexahedra[i_hexahedron][7];
+ }
+ for (size_t i_prism = 0; i_prism < nb_prisms; ++i_prism) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][2];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][3];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][4];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__prisms[i_prism][5];
+ }
+ for (size_t i_pyramid = 0; i_pyramid < nb_pyramids; ++i_pyramid) {
+ cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][0];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][1];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][2];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][3];
+ cell_to_node_list[i_cell_node++] = gmsh_data.__pyramids[i_pyramid][4];
}
- descriptor.cell_type_vector[j] = CellType::Tetrahedron;
- descriptor.cell_number_vector[j] = gmsh_data.__tetrahedra_number[j];
}
- const size_t nb_hexahedra = gmsh_data.__hexahedra.size();
+ for (size_t j = 0; j < nb_tetrahedra; ++j) {
+ cell_type_vector[j] = CellType::Tetrahedron;
+ cell_number_vector[j] = gmsh_data.__tetrahedra_number[j];
+ }
+
for (size_t j = 0; j < nb_hexahedra; ++j) {
const size_t jh = nb_tetrahedra + j;
- descriptor.cell_to_node_vector[jh].resize(8);
- for (int r = 0; r < 8; ++r) {
- descriptor.cell_to_node_vector[jh][r] = gmsh_data.__hexahedra[j][r];
- }
- descriptor.cell_type_vector[jh] = CellType::Hexahedron;
- descriptor.cell_number_vector[jh] = gmsh_data.__hexahedra_number[j];
+
+ cell_type_vector[jh] = CellType::Hexahedron;
+ cell_number_vector[jh] = gmsh_data.__hexahedra_number[j];
}
- const size_t nb_prisms = gmsh_data.__prisms.size();
for (size_t j = 0; j < nb_prisms; ++j) {
const size_t jp = nb_tetrahedra + nb_hexahedra + j;
- descriptor.cell_to_node_vector[jp].resize(6);
- for (int r = 0; r < 6; ++r) {
- descriptor.cell_to_node_vector[jp][r] = gmsh_data.__prisms[j][r];
- }
- descriptor.cell_type_vector[jp] = CellType::Prism;
- descriptor.cell_number_vector[jp] = gmsh_data.__prisms_number[j];
+
+ cell_type_vector[jp] = CellType::Prism;
+ cell_number_vector[jp] = gmsh_data.__prisms_number[j];
}
- const size_t nb_pyramids = gmsh_data.__pyramids.size();
for (size_t j = 0; j < nb_pyramids; ++j) {
const size_t jh = nb_tetrahedra + nb_hexahedra + nb_prisms + j;
- descriptor.cell_to_node_vector[jh].resize(5);
- for (int r = 0; r < 5; ++r) {
- descriptor.cell_to_node_vector[jh][r] = gmsh_data.__pyramids[j][r];
- }
- descriptor.cell_type_vector[jh] = CellType::Pyramid;
- descriptor.cell_number_vector[jh] = gmsh_data.__pyramids_number[j];
+
+ cell_type_vector[jh] = CellType::Pyramid;
+ cell_number_vector[jh] = gmsh_data.__pyramids_number[j];
}
+ descriptor.setCellNumberVector(cell_number_vector);
+ descriptor.setCellTypeVector(cell_type_vector);
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
+
std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
for (unsigned int j = 0; j < gmsh_data.__tetrahedra_ref.size(); ++j) {
const unsigned int elem_number = j;
@@ -429,103 +535,140 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
- const auto& node_number_vector = descriptor.node_number_vector;
+ const auto& node_number_vector = descriptor.nodeNumberVector();
- Array<size_t> face_nb_cell(descriptor.face_number_vector.size());
+ Array<size_t> face_nb_cell(descriptor.faceNumberVector().size());
face_nb_cell.fill(0);
- for (size_t j = 0; j < descriptor.cell_to_face_vector.size(); ++j) {
- for (size_t l = 0; l < descriptor.cell_to_face_vector[j].size(); ++l) {
- face_nb_cell[descriptor.cell_to_face_vector[j][l]] += 1;
+ const auto& cell_to_face_matrix = descriptor.cellToFaceMatrix();
+
+ for (size_t j = 0; j < cell_to_face_matrix.numberOfRows(); ++j) {
+ const auto& cell_face_list = cell_to_face_matrix[j];
+ for (size_t l = 0; l < cell_face_list.size(); ++l) {
+ face_nb_cell[cell_face_list[l]] += 1;
}
}
+ const auto& face_number_vector = descriptor.faceNumberVector();
{
- 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;
+ for (size_t l = 0; l < face_number_vector.size(); ++l) {
+ face_number_id_map[face_number_vector[l]] = l;
}
- Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
+ Assert(face_number_id_map.size() == face_number_vector.size());
return face_number_id_map;
}();
+ const auto& node_to_face_matrix = descriptor.nodeToFaceMatrix();
+ const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
+
+ const auto find_face = [&](std::vector<uint32_t> node_list) {
+ size_t i_node_smallest_number = 0;
+ for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+ if (node_number_vector[node_list[i_node]] < node_number_vector[node_list[i_node_smallest_number]]) {
+ i_node_smallest_number = i_node;
+ }
+ }
+
+ if (i_node_smallest_number != 0) {
+ std::vector<uint64_t> buffer(node_list.size());
+ for (size_t i_node = i_node_smallest_number; i_node < buffer.size(); ++i_node) {
+ buffer[i_node - i_node_smallest_number] = node_list[i_node];
+ }
+ for (size_t i_node = 0; i_node < i_node_smallest_number; ++i_node) {
+ buffer[i_node + node_list.size() - i_node_smallest_number] = node_list[i_node];
+ }
+
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ node_list[i_node] = buffer[i_node];
+ }
+ }
+
+ if (node_number_vector[node_list[1]] > node_number_vector[node_list[node_list.size() - 1]]) {
+ for (size_t i_node = 1; i_node <= (node_list.size() + 1) / 2 - 1; ++i_node) {
+ std::swap(node_list[i_node], node_list[node_list.size() - i_node]);
+ }
+ }
+
+ const auto& face_id_vector = node_to_face_matrix[node_list[0]];
+
+ for (size_t i_face = 0; i_face < face_id_vector.size(); ++i_face) {
+ const FaceId face_id = face_id_vector[i_face];
+ const auto& face_node_id_list = face_to_node_matrix[face_id];
+ if (face_node_id_list.size() == node_list.size()) {
+ bool is_same = true;
+ for (size_t i_node = 1; i_node < face_node_id_list.size(); ++i_node) {
+ is_same &= (face_node_id_list[i_node] == node_list[i_node]);
+ }
+ if (is_same) {
+ return face_id;
+ }
+ }
+ }
+
+ std::stringstream error_msg;
+ error_msg << "face (" << node_list[0];
+ for (size_t i = 1; i < node_list.size(); ++i) {
+ error_msg << ',' << node_list[i];
+ }
+ error_msg << ") not found";
+ throw NormalError(error_msg.str());
+ };
+
+ Array<int> face_number_vector = copy(descriptor.faceNumberVector());
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
for (unsigned int f = 0; f < gmsh_data.__triangles.size(); ++f) {
- const unsigned int face_id = [&] {
- auto i = face_to_id_map.find(
- Face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]},
- node_number_vector));
- if (i == face_to_id_map.end()) {
- throw NormalError("face not found");
- }
- return i->second;
- }();
+ const unsigned int face_id =
+ find_face({gmsh_data.__triangles[f][0], gmsh_data.__triangles[f][1], gmsh_data.__triangles[f][2]});
const unsigned int& ref = gmsh_data.__triangles_ref[f];
ref_faces_map[ref].push_back(face_id);
- if (descriptor.face_number_vector[face_id] != gmsh_data.__triangles_number[f]) {
+ if (face_number_vector[face_id] != gmsh_data.__triangles_number[f]) {
if (auto i_face = face_number_id_map.find(gmsh_data.__triangles_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]);
+ std::swap(face_number_vector[face_id], 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.erase(face_number_vector[face_id]);
+ face_number_id_map.erase(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;
+ face_number_id_map[face_number_vector[face_id]] = face_id;
+ face_number_id_map[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] = gmsh_data.__triangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ face_number_id_map.erase(face_number_vector[face_id]);
+ face_number_vector[face_id] = gmsh_data.__triangles_number[f];
+ face_number_id_map[face_number_vector[face_id]] = face_id;
}
}
}
for (unsigned int f = 0; f < gmsh_data.__quadrangles.size(); ++f) {
- const unsigned int face_id = [&] {
- auto i = face_to_id_map.find(Face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
- gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]},
- node_number_vector));
- if (i == face_to_id_map.end()) {
- throw NormalError("face not found");
- }
- return i->second;
- }();
+ const unsigned int face_id = find_face({gmsh_data.__quadrangles[f][0], gmsh_data.__quadrangles[f][1],
+ gmsh_data.__quadrangles[f][2], gmsh_data.__quadrangles[f][3]});
const unsigned int& ref = gmsh_data.__quadrangles_ref[f];
ref_faces_map[ref].push_back(face_id);
-
- if (descriptor.face_number_vector[face_id] != gmsh_data.__quadrangles_number[f]) {
+ if (face_number_vector[face_id] != gmsh_data.__quadrangles_number[f]) {
if (auto i_face = face_number_id_map.find(gmsh_data.__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]);
+ std::swap(face_number_vector[face_id], 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.erase(face_number_vector[face_id]);
+ face_number_id_map.erase(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;
+ face_number_id_map[face_number_vector[face_id]] = face_id;
+ face_number_id_map[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] = gmsh_data.__quadrangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
+ face_number_id_map.erase(face_number_vector[face_id]);
+ face_number_vector[face_id] = gmsh_data.__quadrangles_number[f];
+ face_number_id_map[face_number_vector[face_id]] = face_id;
}
}
}
+ descriptor.setFaceNumberVector(face_number_vector);
for (const auto& ref_face_list : ref_faces_map) {
Array<FaceId> face_list(ref_face_list.second.size());
@@ -556,68 +699,76 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(descriptor);
{
- Array<bool> is_boundary_edge(descriptor.edge_number_vector.size());
+ Array<bool> is_boundary_edge(descriptor.edgeNumberVector().size());
is_boundary_edge.fill(false);
+ const auto& face_to_edge_matrix = descriptor.faceToEdgeMatrix();
+
for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
if (face_nb_cell[i_face] == 1) {
- for (size_t node_id : descriptor.face_to_edge_vector[i_face]) {
- is_boundary_edge[node_id] = true;
+ auto face_edge_list = face_to_edge_matrix[i_face];
+ for (size_t i_edge = 0; i_edge < face_edge_list.size(); ++i_edge) {
+ is_boundary_edge[face_edge_list[i_edge]] = true;
}
}
}
- 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;
+ const auto& node_to_edge_matrix = descriptor.nodeToEdgeMatrix();
+ const auto& edge_to_node_matrix = descriptor.edgeToNodeMatrix();
+
+ const auto find_edge = [&](uint32_t node0, uint32_t node1) {
+ if (node_number_vector[node0] > node_number_vector[node1]) {
+ std::swap(node0, node1);
}
- return edge_to_id_map;
- }();
+ const auto& edge_id_vector = node_to_edge_matrix[node0];
+
+ for (size_t i_edge = 0; i_edge < edge_id_vector.size(); ++i_edge) {
+ const EdgeId edge_id = edge_id_vector[i_edge];
+ if (edge_to_node_matrix[edge_id][1] == node1) {
+ return edge_id;
+ }
+ }
+
+ std::stringstream error_msg;
+ error_msg << "edge (" << node0 << ',' << node1 << ") not found";
+ throw NormalError(error_msg.str());
+ };
std::unordered_map<int, EdgeId> edge_number_id_map = [&] {
+ const auto& edge_number_vector = descriptor.edgeNumberVector();
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;
+ for (size_t l = 0; l < edge_number_vector.size(); ++l) {
+ edge_number_id_map[edge_number_vector[l]] = l;
}
- Assert(edge_number_id_map.size() == descriptor.edge_number_vector.size());
+ Assert(edge_number_id_map.size() == edge_number_vector.size());
return edge_number_id_map;
}();
+ Array<int> edge_number_vector = copy(descriptor.edgeNumberVector());
std::map<unsigned int, std::vector<unsigned int>> ref_edges_map;
for (unsigned int e = 0; e < gmsh_data.__edges.size(); ++e) {
- const unsigned int edge_id = [&] {
- auto i = edge_to_id_map.find(Edge({gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]}, node_number_vector));
- if (i == edge_to_id_map.end()) {
- std::stringstream error_msg;
- error_msg << "edge " << gmsh_data.__edges[e][0] << " not found";
- throw NormalError(error_msg.str());
- }
- return i->second;
- }();
- const unsigned int& ref = gmsh_data.__edges_ref[e];
+ const unsigned int edge_id = find_edge(gmsh_data.__edges[e][0], gmsh_data.__edges[e][1]);
+ const unsigned int& ref = gmsh_data.__edges_ref[e];
ref_edges_map[ref].push_back(edge_id);
- if (descriptor.edge_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
+ if (edge_number_vector[edge_id] != gmsh_data.__edges_number[e]) {
if (auto i_edge = edge_number_id_map.find(gmsh_data.__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]);
+ std::swap(edge_number_vector[edge_id], 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.erase(edge_number_vector[edge_id]);
+ edge_number_id_map.erase(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;
+ edge_number_id_map[edge_number_vector[edge_id]] = edge_id;
+ edge_number_id_map[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] = gmsh_data.__edges_number[e];
- edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_id;
+ edge_number_id_map.erase(edge_number_vector[edge_id]);
+ edge_number_vector[edge_id] = gmsh_data.__edges_number[e];
+ edge_number_id_map[edge_number_vector[edge_id]] = edge_id;
}
}
}
+ descriptor.setEdgeNumberVector(edge_number_vector);
for (const auto& ref_edge_list : ref_edges_map) {
Array<EdgeId> edge_list(ref_edge_list.second.size());
@@ -645,12 +796,16 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
}
}
- Array<bool> is_boundary_node(descriptor.node_number_vector.size());
+ Array<bool> is_boundary_node(node_number_vector.size());
is_boundary_node.fill(false);
+ const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
+
for (size_t i_face = 0; i_face < face_nb_cell.size(); ++i_face) {
if (face_nb_cell[i_face] == 1) {
- for (size_t node_id : descriptor.face_to_node_vector[i_face]) {
+ const auto& face_node_list = face_to_node_matrix[i_face];
+ for (size_t i_node = 0; i_node < face_node_list.size(); ++i_node) {
+ const NodeId node_id = face_node_list[i_node];
is_boundary_node[node_id] = true;
}
}
@@ -688,17 +843,29 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
descriptor.addRefItemList(RefNodeList(ref_id, point_list, is_boundary));
}
- 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());
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(nb_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
+
+ descriptor.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector(descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::rank());
+ return face_owner_vector;
+ }());
+
+ descriptor.setEdgeOwnerVector([&] {
+ Array<int> edge_owner_vector(descriptor.edgeNumberVector().size());
+ edge_owner_vector.fill(parallel::rank());
+ return edge_owner_vector;
+ }());
+
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(node_number_vector.size());
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
m_connectivity = Connectivity3D::build(descriptor);
}
@@ -958,255 +1125,6 @@ GmshReader::__readPeriodic2_2()
}
}
-// std::shared_ptr<IConnectivity>
-// GmshReader::_buildConnectivity3D(const size_t nb_cells)
-// {
-// ConnectivityDescriptor descriptor;
-
-// descriptor.node_number_vector = m_mesh_data.__verticesNumbers;
-// descriptor.cell_type_vector.resize(nb_cells);
-// descriptor.cell_number_vector.resize(nb_cells);
-// descriptor.cell_to_node_vector.resize(nb_cells);
-
-// const size_t nb_tetrahedra = m_mesh_data.__tetrahedra.size();
-// for (size_t j = 0; j < nb_tetrahedra; ++j) {
-// descriptor.cell_to_node_vector[j].resize(4);
-// for (int r = 0; r < 4; ++r) {
-// descriptor.cell_to_node_vector[j][r] = m_mesh_data.__tetrahedra[j][r];
-// }
-// descriptor.cell_type_vector[j] = CellType::Tetrahedron;
-// descriptor.cell_number_vector[j] = m_mesh_data.__tetrahedra_number[j];
-// }
-// const size_t nb_hexahedra = m_mesh_data.__hexahedra.size();
-// for (size_t j = 0; j < nb_hexahedra; ++j) {
-// const size_t jh = nb_tetrahedra + j;
-// descriptor.cell_to_node_vector[jh].resize(8);
-// for (int r = 0; r < 8; ++r) {
-// descriptor.cell_to_node_vector[jh][r] = m_mesh_data.__hexahedra[j][r];
-// }
-// descriptor.cell_type_vector[jh] = CellType::Hexahedron;
-// descriptor.cell_number_vector[jh] = m_mesh_data.__hexahedra_number[j];
-// }
-
-// std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
-// for (unsigned int r = 0; r < m_mesh_data.__tetrahedra_ref.size(); ++r) {
-// const unsigned int elem_number = m_mesh_data.__tetrahedra_ref[r];
-// const unsigned int& ref = m_mesh_data.__tetrahedra_ref[r];
-// ref_cells_map[ref].push_back(elem_number);
-// }
-
-// for (unsigned int j = 0; j < m_mesh_data.__hexahedra_ref.size(); ++j) {
-// const size_t elem_number = nb_tetrahedra + j;
-// const unsigned int& ref = m_mesh_data.__hexahedra_ref[j];
-// ref_cells_map[ref].push_back(elem_number);
-// }
-
-// for (const auto& ref_cell_list : ref_cells_map) {
-// Array<CellId> cell_list(ref_cell_list.second.size());
-// for (size_t j = 0; j < ref_cell_list.second.size(); ++j) {
-// cell_list[j] = ref_cell_list.second[j];
-// }
-// const PhysicalRefId& physical_ref_id = m_mesh_data.m_physical_ref_map.at(ref_cell_list.first);
-// descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
-// }
-
-// ConnectivityBuilderBase::_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 < m_mesh_data.__triangles.size(); ++f) {
-// const unsigned int face_id = [&] {
-// auto i = face_to_id_map.find(
-// Face({m_mesh_data.__triangles[f][0], m_mesh_data.__triangles[f][1], m_mesh_data.__triangles[f][2]},
-// node_number_vector));
-// if (i == face_to_id_map.end()) {
-// throw NormalError("face not found");
-// }
-// return i->second;
-// }();
-
-// const unsigned int& ref = m_mesh_data.__triangles_ref[f];
-// ref_faces_map[ref].push_back(face_id);
-
-// if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
-// if (auto i_face = face_number_id_map.find(m_mesh_data.__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] = m_mesh_data.__quadrangles_number[f];
-// face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
-// }
-// }
-// }
-
-// for (unsigned int f = 0; f < m_mesh_data.__quadrangles.size(); ++f) {
-// const unsigned int face_id = [&] {
-// auto i = face_to_id_map.find(Face({m_mesh_data.__quadrangles[f][0], m_mesh_data.__quadrangles[f][1],
-// m_mesh_data.__quadrangles[f][2], m_mesh_data.__quadrangles[f][3]},
-// node_number_vector));
-// if (i == face_to_id_map.end()) {
-// throw NormalError("face not found");
-// }
-// return i->second;
-// }();
-
-// const unsigned int& ref = m_mesh_data.__quadrangles_ref[f];
-// ref_faces_map[ref].push_back(face_id);
-
-// if (descriptor.face_number_vector[face_id] != m_mesh_data.__quadrangles_number[f]) {
-// if (auto i_face = face_number_id_map.find(m_mesh_data.__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] = m_mesh_data.__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_mesh_data.m_physical_ref_map.at(ref_face_list.first);
-// descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
-// }
-// }
-
-// ConnectivityBuilderBase::_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 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 < m_mesh_data.__edges.size(); ++e) {
-// const unsigned int edge_id = [&] {
-// auto i = edge_to_id_map.find(Edge({m_mesh_data.__edges[e][0], m_mesh_data.__edges[e][1]},
-// node_number_vector)); if (i == edge_to_id_map.end()) {
-// std::stringstream error_msg;
-// error_msg << "edge " << m_mesh_data.__edges[e][0] << " not found";
-// throw NormalError(error_msg.str());
-// }
-// return i->second;
-// }();
-// const unsigned int& ref = m_mesh_data.__edges_ref[e];
-// ref_edges_map[ref].push_back(edge_id);
-
-// if (descriptor.edge_number_vector[edge_id] != m_mesh_data.__edges_number[e]) {
-// if (auto i_edge = edge_number_id_map.find(m_mesh_data.__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] = m_mesh_data.__edges_number[e];
-// edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_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_mesh_data.m_physical_ref_map.at(ref_edge_list.first);
-// descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
-// }
-// }
-
-// std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
-// for (unsigned int r = 0; r < m_mesh_data.__points.size(); ++r) {
-// const unsigned int point_number = m_mesh_data.__points[r];
-// const unsigned int& ref = m_mesh_data.__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_mesh_data.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());
-
-// return Connectivity3D::build(descriptor);
-// }
-
void
GmshReader::__proceedData()
{
diff --git a/src/mesh/LogicalConnectivityBuilder.cpp b/src/mesh/LogicalConnectivityBuilder.cpp
index 41092faaf90a9a3f9a24ac1c561304373ba748d0..dbc41d0681e9ceda822a5896af97ace2e3572f03 100644
--- a/src/mesh/LogicalConnectivityBuilder.cpp
+++ b/src/mesh/LogicalConnectivityBuilder.cpp
@@ -141,28 +141,27 @@ void
LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
ConnectivityDescriptor& 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;
+ const auto& node_number_vector = descriptor.nodeNumberVector();
+ const auto& node_to_edge_matrix = descriptor.nodeToEdgeMatrix();
+ const auto& edge_to_node_matrix = descriptor.edgeToNodeMatrix();
+
+ const auto find_edge = [&](uint32_t node0, uint32_t node1) {
+ if (node_number_vector[node0] > node_number_vector[node1]) {
+ std::swap(node0, node1);
}
- return edge_to_id_map;
- }();
+ const auto& node_edge_list = node_to_edge_matrix[node0];
- 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;
+ for (size_t i_edge = 0; i_edge < node_edge_list.size(); ++i_edge) {
+ const EdgeId edge_id = node_edge_list[i_edge];
+ if (edge_to_node_matrix[edge_id][1] == node1) {
+ return edge_id;
+ }
}
- Assert(edge_number_id_map.size() == descriptor.edge_number_vector.size());
- return edge_number_id_map;
- }();
+ throw UnexpectedError("Cannot find edge");
+ };
const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
- const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+ const auto get_node_id = [&](const TinyVector<3, uint64_t>& node_logic_id) {
return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
};
@@ -172,13 +171,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
Array<EdgeId> boundary_edges(cell_size[2]);
size_t l = 0;
for (size_t k = 0; k < cell_size[2]; ++k) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
- auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
- Assert(i_edge != edge_to_id_map.end());
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
- boundary_edges[l++] = i_edge->second;
+ boundary_edges[l++] = find_edge(node_0_id, node_1_id);
}
Assert(l == cell_size[2]);
descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -196,13 +192,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
Array<EdgeId> boundary_edges(cell_size[1]);
size_t l = 0;
for (size_t j = 0; j < cell_size[1]; ++j) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
-
- auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
- Assert(i_edge != edge_to_id_map.end());
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
- boundary_edges[l++] = i_edge->second;
+ boundary_edges[l++] = find_edge(node_0_id, node_1_id);
}
Assert(l == cell_size[1]);
descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -220,13 +213,10 @@ LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>
Array<EdgeId> boundary_edges(cell_size[0]);
size_t l = 0;
for (size_t i = 0; i < cell_size[0]; ++i) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
- auto i_edge = edge_to_id_map.find(Edge{{node_0_id, node_1_id}, descriptor.node_number_vector});
- Assert(i_edge != edge_to_id_map.end());
-
- boundary_edges[l++] = i_edge->second;
+ boundary_edges[l++] = find_edge(node_0_id, node_1_id);
}
Assert(l == cell_size[0]);
descriptor.addRefItemList(RefEdgeList{RefId{ref_id, ref_name}, boundary_edges, true});
@@ -256,6 +246,8 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
return cell_logic_id[0] * cell_size[1] + cell_logic_id[1];
};
+ const auto& cell_to_face_matrix = descriptor.cellToFaceMatrix();
+
{ // xmin
const size_t i = 0;
Array<FaceId> boundary_faces(cell_size[1]);
@@ -263,7 +255,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
constexpr size_t left_face = 3;
const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
- const size_t face_id = descriptor.cell_to_face_vector[cell_id][left_face];
+ const size_t face_id = cell_to_face_matrix[cell_id][left_face];
boundary_faces[j] = face_id;
}
@@ -277,7 +269,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
constexpr size_t right_face = 1;
const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
- const size_t face_id = descriptor.cell_to_face_vector[cell_id][right_face];
+ const size_t face_id = cell_to_face_matrix[cell_id][right_face];
boundary_faces[j] = face_id;
}
@@ -291,7 +283,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
constexpr size_t bottom_face = 0;
const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
- const size_t face_id = descriptor.cell_to_face_vector[cell_id][bottom_face];
+ const size_t face_id = cell_to_face_matrix[cell_id][bottom_face];
boundary_faces[i] = face_id;
}
@@ -305,7 +297,7 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(
constexpr size_t top_face = 2;
const size_t cell_id = cell_number(TinyVector<2, uint64_t>{i, j});
- const size_t face_id = descriptor.cell_to_face_vector[cell_id][top_face];
+ const size_t face_id = cell_to_face_matrix[cell_id][top_face];
boundary_faces[i] = face_id;
}
@@ -318,28 +310,53 @@ void
LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
ConnectivityDescriptor& descriptor)
{
- using Face = ConnectivityFace<3>;
+ const auto& node_number_vector = descriptor.nodeNumberVector();
+ const auto& face_to_node_matrix = descriptor.faceToNodeMatrix();
+ const auto& node_to_face_matrix = descriptor.nodeToFaceMatrix();
+
+ const auto find_face = [&](std::array<uint32_t, 4> node_list) {
+ size_t i_node_smallest_number = 0;
+ for (size_t i_node = 1; i_node < node_list.size(); ++i_node) {
+ if (node_number_vector[node_list[i_node]] < node_number_vector[node_list[i_node_smallest_number]]) {
+ i_node_smallest_number = i_node;
+ }
+ }
+
+ if (i_node_smallest_number != 0) {
+ std::array<uint64_t, 4> buffer;
+ for (size_t i_node = i_node_smallest_number; i_node < buffer.size(); ++i_node) {
+ buffer[i_node - i_node_smallest_number] = node_list[i_node];
+ }
+ for (size_t i_node = 0; i_node < i_node_smallest_number; ++i_node) {
+ buffer[i_node + node_list.size() - i_node_smallest_number] = node_list[i_node];
+ }
- 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, descriptor.node_number_vector)] = l;
+ for (size_t i_node = 0; i_node < node_list.size(); ++i_node) {
+ node_list[i_node] = buffer[i_node];
+ }
}
- return face_to_id_map;
- }();
- const 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;
+ if (node_number_vector[node_list[1]] > node_number_vector[node_list[node_list.size() - 1]]) {
+ for (size_t i_node = 1; i_node <= (node_list.size() + 1) / 2 - 1; ++i_node) {
+ std::swap(node_list[i_node], node_list[node_list.size() - i_node]);
+ }
}
- Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
- return face_number_id_map;
- }();
+
+ const auto& node_face_list = node_to_face_matrix[node_list[0]];
+
+ for (size_t i_face = 0; i_face < node_face_list.size(); ++i_face) {
+ const FaceId face_id = node_face_list[i_face];
+ const auto& face_node_list = face_to_node_matrix[face_id];
+ if ((face_node_list[1] == node_list[1]) and (face_node_list[2] == node_list[2]) and
+ (face_node_list[3] == node_list[3])) {
+ return face_id;
+ }
+ }
+ throw UnexpectedError("Cannot find edge");
+ };
const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
- const auto node_number = [&](const TinyVector<3, uint64_t>& node_logic_id) {
+ const auto get_node_id = [&](const TinyVector<3, uint64_t>& node_logic_id) {
return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
};
@@ -349,16 +366,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
size_t l = 0;
for (size_t j = 0; j < cell_size[1]; ++j) {
for (size_t k = 0; k < cell_size[2]; ++k) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
- const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k + 1});
- const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
+ const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k + 1});
+ const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
- auto i_face =
- face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
- Assert(i_face != face_to_id_map.end());
-
- boundary_faces[l++] = i_face->second;
+ boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
}
}
Assert(l == cell_size[1] * cell_size[2]);
@@ -375,16 +388,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
size_t l = 0;
for (size_t i = 0; i < cell_size[0]; ++i) {
for (size_t k = 0; k < cell_size[2]; ++k) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
- const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k + 1});
- const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j, k + 1});
-
- auto i_face =
- face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
- Assert(i_face != face_to_id_map.end());
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k + 1});
+ const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j, k + 1});
- boundary_faces[l++] = i_face->second;
+ boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
}
}
Assert(l == cell_size[0] * cell_size[2]);
@@ -401,16 +410,12 @@ LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>
size_t l = 0;
for (size_t i = 0; i < cell_size[0]; ++i) {
for (size_t j = 0; j < cell_size[1]; ++j) {
- const uint32_t node_0_id = node_number(TinyVector<3, uint64_t>{i, j, k});
- const uint32_t node_1_id = node_number(TinyVector<3, uint64_t>{i + 1, j, k});
- const uint32_t node_2_id = node_number(TinyVector<3, uint64_t>{i + 1, j + 1, k});
- const uint32_t node_3_id = node_number(TinyVector<3, uint64_t>{i, j + 1, k});
+ const uint32_t node_0_id = get_node_id(TinyVector<3, uint64_t>{i, j, k});
+ const uint32_t node_1_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j, k});
+ const uint32_t node_2_id = get_node_id(TinyVector<3, uint64_t>{i + 1, j + 1, k});
+ const uint32_t node_3_id = get_node_id(TinyVector<3, uint64_t>{i, j + 1, k});
- auto i_face =
- face_to_id_map.find(Face{{node_0_id, node_1_id, node_2_id, node_3_id}, descriptor.node_number_vector});
- Assert(i_face != face_to_id_map.end());
-
- boundary_faces[l++] = i_face->second;
+ boundary_faces[l++] = find_face({node_0_id, node_1_id, node_2_id, node_3_id});
}
}
Assert(l == cell_size[0] * cell_size[1]);
@@ -431,36 +436,50 @@ LogicalConnectivityBuilder::_buildConnectivity(
const size_t number_of_nodes = cell_size[0] + 1;
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector.resize(number_of_nodes);
- for (size_t i = 0; i < number_of_nodes; ++i) {
- descriptor.node_number_vector[i] = i;
- }
+ descriptor.setNodeNumberVector([&] {
+ Array<int> node_number_vector(number_of_nodes);
+ parallel_for(
+ number_of_nodes, PUGS_LAMBDA(const size_t i) { node_number_vector[i] = i; });
+ return node_number_vector;
+ }());
+
+ descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector(number_of_cells);
+ for (size_t i = 0; i < number_of_cells; ++i) {
+ cell_number_vector[i] = i;
+ }
+ return cell_number_vector;
+ }());
- descriptor.cell_number_vector.resize(number_of_cells);
+ Array<CellType> cell_type_vector(number_of_cells);
+ cell_type_vector.fill(CellType::Line);
+ descriptor.setCellTypeVector(cell_type_vector);
+
+ Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+ for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+ cell_to_node_row_map[i] = 2 * i;
+ }
+ Array<unsigned int> cell_to_node_list(2 * number_of_cells);
for (size_t i = 0; i < number_of_cells; ++i) {
- descriptor.cell_number_vector[i] = i;
+ cell_to_node_list[2 * i] = i;
+ cell_to_node_list[2 * i + 1] = i + 1;
}
- descriptor.cell_type_vector.resize(number_of_cells);
- std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Line);
-
- descriptor.cell_to_node_vector.resize(number_of_cells);
- constexpr size_t nb_node_per_cell = 2;
- for (size_t j = 0; j < number_of_cells; ++j) {
- descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
- for (size_t r = 0; r < nb_node_per_cell; ++r) {
- descriptor.cell_to_node_vector[j][0] = j;
- descriptor.cell_to_node_vector[j][1] = j + 1;
- }
- }
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list));
this->_buildBoundaryNodeList(cell_size, descriptor);
- descriptor.cell_owner_vector.resize(number_of_cells);
- std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(number_of_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(number_of_nodes);
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
m_connectivity = Connectivity1D::build(descriptor);
}
@@ -478,18 +497,23 @@ LogicalConnectivityBuilder::_buildConnectivity(
const size_t number_of_nodes = node_size[0] * node_size[1];
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector.resize(number_of_nodes);
- for (size_t i = 0; i < number_of_nodes; ++i) {
- descriptor.node_number_vector[i] = i;
- }
-
- descriptor.cell_number_vector.resize(number_of_cells);
- for (size_t i = 0; i < number_of_cells; ++i) {
- descriptor.cell_number_vector[i] = i;
- }
-
- descriptor.cell_type_vector.resize(number_of_cells);
- std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Quadrangle);
+ descriptor.setNodeNumberVector([&] {
+ Array<int> node_number_vector(number_of_nodes);
+ parallel_for(
+ number_of_nodes, PUGS_LAMBDA(const size_t i) { node_number_vector[i] = i; });
+ return node_number_vector;
+ }());
+
+ descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector(number_of_cells);
+ parallel_for(
+ number_of_cells, PUGS_LAMBDA(size_t i) { cell_number_vector[i] = i; });
+ return cell_number_vector;
+ }());
+
+ Array<CellType> cell_type_vector(number_of_cells);
+ cell_type_vector.fill(CellType::Quadrangle);
+ descriptor.setCellTypeVector(cell_type_vector);
const auto node_number = [&](const TinyVector<Dimension, uint64_t> node_logic_id) {
return node_logic_id[0] * node_size[1] + node_logic_id[1];
@@ -501,32 +525,44 @@ LogicalConnectivityBuilder::_buildConnectivity(
return TinyVector<Dimension, uint64_t>{j0, j1};
};
- descriptor.cell_to_node_vector.resize(number_of_cells);
- constexpr size_t nb_node_per_cell = 1 << Dimension;
+ Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+ for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+ cell_to_node_row_map[i] = 4 * i;
+ }
+ Array<unsigned int> cell_to_node_list(4 * number_of_cells);
for (size_t j = 0; j < number_of_cells; ++j) {
TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
- descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
- for (size_t r = 0; r < nb_node_per_cell; ++r) {
- descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
- descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
- descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
- descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
- }
+
+ cell_to_node_list[4 * j] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0});
+ cell_to_node_list[4 * j + 1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0});
+ cell_to_node_list[4 * j + 2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1});
+ cell_to_node_list[4 * j + 3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1});
}
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list));
+
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
this->_buildBoundaryNodeList(cell_size, descriptor);
this->_buildBoundaryFaceList(cell_size, descriptor);
- descriptor.cell_owner_vector.resize(number_of_cells);
- std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(number_of_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
- 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.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector(descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::rank());
+ return face_owner_vector;
+ }());
- descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(descriptor.nodeNumberVector().size());
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
m_connectivity = Connectivity<Dimension>::build(descriptor);
}
@@ -560,18 +596,23 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
const size_t number_of_nodes = count_items(node_size);
ConnectivityDescriptor descriptor;
- descriptor.node_number_vector.resize(number_of_nodes);
- for (size_t i = 0; i < number_of_nodes; ++i) {
- descriptor.node_number_vector[i] = i;
- }
-
- descriptor.cell_number_vector.resize(number_of_cells);
- for (size_t i = 0; i < number_of_cells; ++i) {
- descriptor.cell_number_vector[i] = i;
- }
-
- descriptor.cell_type_vector.resize(number_of_cells);
- std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Hexahedron);
+ descriptor.setNodeNumberVector([&] {
+ Array<int> node_number_vector(number_of_nodes);
+ parallel_for(
+ number_of_nodes, PUGS_LAMBDA(const size_t i) { node_number_vector[i] = i; });
+ return node_number_vector;
+ }());
+
+ descriptor.setCellNumberVector([&] {
+ Array<int> cell_number_vector(number_of_cells);
+ parallel_for(
+ number_of_cells, PUGS_LAMBDA(size_t i) { cell_number_vector[i] = i; });
+ return cell_number_vector;
+ }());
+
+ Array<CellType> cell_type_vector(number_of_cells);
+ cell_type_vector.fill(CellType::Hexahedron);
+ descriptor.setCellTypeVector(cell_type_vector);
const auto cell_logic_id = [&](size_t j) {
const size_t slice1 = cell_size[1] * cell_size[2];
@@ -586,24 +627,26 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
};
- descriptor.cell_to_node_vector.resize(number_of_cells);
- constexpr size_t nb_node_per_cell = 1 << Dimension;
+ Array<unsigned int> cell_to_node_row_map(number_of_cells + 1);
+ for (size_t i = 0; i < cell_to_node_row_map.size(); ++i) {
+ cell_to_node_row_map[i] = 8 * i;
+ }
+ Array<unsigned int> cell_to_node_list(8 * number_of_cells);
for (size_t j = 0; j < number_of_cells; ++j) {
TinyVector<Dimension, size_t> cell_index = cell_logic_id(j);
- descriptor.cell_to_node_vector[j].resize(nb_node_per_cell);
- for (size_t r = 0; r < nb_node_per_cell; ++r) {
- static_assert(Dimension == 3, "unexpected dimension");
- descriptor.cell_to_node_vector[j][0] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
- descriptor.cell_to_node_vector[j][1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
- descriptor.cell_to_node_vector[j][2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
- descriptor.cell_to_node_vector[j][3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
- descriptor.cell_to_node_vector[j][4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
- descriptor.cell_to_node_vector[j][5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
- descriptor.cell_to_node_vector[j][6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
- descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
- }
+
+ cell_to_node_list[8 * j] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 0});
+ cell_to_node_list[8 * j + 1] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 0});
+ cell_to_node_list[8 * j + 2] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 0});
+ cell_to_node_list[8 * j + 3] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 0});
+ cell_to_node_list[8 * j + 4] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 0, 1});
+ cell_to_node_list[8 * j + 5] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 0, 1});
+ cell_to_node_list[8 * j + 6] = node_number(cell_index + TinyVector<Dimension, uint64_t>{1, 1, 1});
+ cell_to_node_list[8 * j + 7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
}
+ descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row_map, cell_to_node_list));
+
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
@@ -611,17 +654,29 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
this->_buildBoundaryEdgeList(cell_size, descriptor);
this->_buildBoundaryFaceList(cell_size, descriptor);
- descriptor.cell_owner_vector.resize(number_of_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());
+ descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(number_of_cells);
+ cell_owner_vector.fill(parallel::rank());
+ return cell_owner_vector;
+ }());
+
+ descriptor.setFaceOwnerVector([&] {
+ Array<int> face_owner_vector(descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::rank());
+ return face_owner_vector;
+ }());
+
+ descriptor.setEdgeOwnerVector([&] {
+ Array<int> edge_owner_vector(descriptor.edgeNumberVector().size());
+ edge_owner_vector.fill(parallel::rank());
+ return edge_owner_vector;
+ }());
+
+ descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(descriptor.nodeNumberVector().size());
+ node_owner_vector.fill(parallel::rank());
+ return node_owner_vector;
+ }());
m_connectivity = Connectivity<Dimension>::build(descriptor);
}
diff --git a/src/mesh/MedianDualConnectivityBuilder.cpp b/src/mesh/MedianDualConnectivityBuilder.cpp
index c6770d05ee1d93891d560bd1b1373fd2e68da5eb..1986c4a826a0a3ba2754e6204dafb4f9151bd5f7 100644
--- a/src/mesh/MedianDualConnectivityBuilder.cpp
+++ b/src/mesh/MedianDualConnectivityBuilder.cpp
@@ -92,33 +92,37 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
Assert(i_boundary_node == primal_number_of_boundary_nodes);
}
- dual_descriptor.node_number_vector.resize(dual_number_of_nodes);
+ Array<int> node_number_vector(dual_number_of_nodes);
{
parallel_for(m_primal_cell_to_dual_node_map.size(), [&](size_t i) {
- const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
- dual_descriptor.node_number_vector[dual_node_id] = primal_cell_number[primal_cell_id];
+ const auto [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i];
+ node_number_vector[dual_node_id] = primal_cell_number[primal_cell_id];
});
const size_t face_number_shift = max(primal_cell_number) + 1;
parallel_for(primal_number_of_faces, [&](size_t i) {
- const auto [primal_face_id, dual_node_id] = m_primal_face_to_dual_node_map[i];
- dual_descriptor.node_number_vector[dual_node_id] = primal_face_number[primal_face_id] + face_number_shift;
+ const auto [primal_face_id, dual_node_id] = m_primal_face_to_dual_node_map[i];
+ node_number_vector[dual_node_id] = primal_face_number[primal_face_id] + face_number_shift;
});
const size_t node_number_shift = face_number_shift + max(primal_face_number) + 1;
parallel_for(m_primal_boundary_node_to_dual_node_map.size(), [&](size_t i) {
- const auto [primal_node_id, dual_node_id] = m_primal_boundary_node_to_dual_node_map[i];
- dual_descriptor.node_number_vector[dual_node_id] = primal_node_number[primal_node_id] + node_number_shift;
+ const auto [primal_node_id, dual_node_id] = m_primal_boundary_node_to_dual_node_map[i];
+ node_number_vector[dual_node_id] = primal_node_number[primal_node_id] + node_number_shift;
});
}
+ dual_descriptor.setNodeNumberVector(node_number_vector);
- dual_descriptor.cell_number_vector.resize(dual_number_of_cells);
- parallel_for(dual_number_of_cells, [&](size_t i) {
- const auto [primal_node_id, dual_cell_id] = m_primal_node_to_dual_cell_map[i];
- dual_descriptor.cell_number_vector[dual_cell_id] = primal_node_number[primal_node_id];
- });
+ {
+ Array<int> cell_number_vector(dual_number_of_cells);
+ parallel_for(dual_number_of_cells, [&](size_t i) {
+ const auto [primal_node_id, dual_cell_id] = m_primal_node_to_dual_cell_map[i];
+ cell_number_vector[dual_cell_id] = primal_node_number[primal_node_id];
+ });
+ dual_descriptor.setCellNumberVector(cell_number_vector);
+ }
- dual_descriptor.cell_type_vector.resize(dual_number_of_cells);
+ Array<CellType> cell_type_vector(dual_number_of_cells);
const auto& primal_node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
@@ -127,13 +131,13 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
const auto& primal_node_cell_list = primal_node_to_cell_matrix[node_id];
if (primal_node_cell_list.size() == 1) {
- dual_descriptor.cell_type_vector[i_dual_cell] = CellType::Quadrangle;
+ cell_type_vector[i_dual_cell] = CellType::Quadrangle;
} else {
- dual_descriptor.cell_type_vector[i_dual_cell] = CellType::Polygon;
+ cell_type_vector[i_dual_cell] = CellType::Polygon;
}
});
+ dual_descriptor.setCellTypeVector(cell_type_vector);
- dual_descriptor.cell_to_node_vector.resize(dual_number_of_cells);
const auto& primal_cell_to_face_matrix = primal_connectivity.cellToFaceMatrix();
const auto& primal_node_to_face_matrix = primal_connectivity.nodeToFaceMatrix();
const auto& primal_face_to_cell_matrix = primal_connectivity.faceToCellMatrix();
@@ -170,16 +174,35 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
// LCOV_EXCL_STOP
};
+ Array<unsigned int> cell_to_node_row(dual_number_of_cells + 1);
+ cell_to_node_row[0] = 0;
+ {
+ for (NodeId node_id = 0; node_id < primal_number_of_nodes; ++node_id) {
+ // const size_t i_dual_cell = node_id;
+ const auto& primal_node_to_cell_list = primal_node_to_cell_matrix[node_id];
+ const auto& primal_node_to_face_list = primal_node_to_face_matrix[node_id];
+
+ if (primal_node_to_cell_list.size() != primal_node_to_face_list.size()) {
+ // boundary cell
+ cell_to_node_row[node_id + 1] =
+ cell_to_node_row[node_id] + 1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size();
+ } else {
+ // inner cell
+ cell_to_node_row[node_id + 1] =
+ cell_to_node_row[node_id] + primal_node_to_cell_list.size() + primal_node_to_face_list.size();
+ }
+ }
+ }
+
+ Array<unsigned int> cell_to_node_list(cell_to_node_row[cell_to_node_row.size() - 1]);
parallel_for(primal_number_of_nodes, [&](NodeId node_id) {
- const size_t i_dual_cell = node_id;
const auto& primal_node_to_cell_list = primal_node_to_cell_matrix[node_id];
const auto& primal_node_to_face_list = primal_node_to_face_matrix[node_id];
- auto& dual_cell_node_list = dual_descriptor.cell_to_node_vector[i_dual_cell];
+ size_t i_dual_cell_node = cell_to_node_row[node_id];
if (primal_node_to_cell_list.size() != primal_node_to_face_list.size()) {
// boundary cell
- dual_cell_node_list.reserve(1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size());
auto [face_id, cell_id] = [&]() -> std::pair<FaceId, CellId> {
for (size_t i_face = 0; i_face < primal_node_to_face_list.size(); ++i_face) {
@@ -200,24 +223,24 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
// LCOV_EXCL_STOP
}();
- dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
- dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+ cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+ cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
face_id = next_face(cell_id, face_id, node_id);
while (primal_face_to_cell_matrix[face_id].size() > 1) {
- dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
- cell_id = next_cell(cell_id, face_id);
- dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+ cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+
+ cell_id = next_cell(cell_id, face_id);
+ cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
+
face_id = next_face(cell_id, face_id, node_id);
}
- dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
- dual_cell_node_list.push_back(node_to_dual_node_correpondance[node_id]);
+ cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+ cell_to_node_list[i_dual_cell_node++] = node_to_dual_node_correpondance[node_id];
- Assert(dual_cell_node_list.size() == 1 + primal_node_to_cell_list.size() + primal_node_to_face_list.size());
} else {
// inner cell
- dual_cell_node_list.reserve(primal_node_to_cell_list.size() + primal_node_to_face_list.size());
auto [face_id, cell_id] = [&]() -> std::pair<FaceId, CellId> {
const FaceId face_id = primal_node_to_face_list[0];
@@ -237,14 +260,16 @@ MedianDualConnectivityBuilder::_buildConnectivityDescriptor<2>(const Connectivit
const FaceId first_face_id = face_id;
do {
- dual_cell_node_list.push_back(m_primal_face_to_dual_node_map[face_id].second);
- dual_cell_node_list.push_back(m_primal_cell_to_dual_node_map[cell_id].second);
+ cell_to_node_list[i_dual_cell_node++] = m_primal_face_to_dual_node_map[face_id].second;
+ cell_to_node_list[i_dual_cell_node++] = m_primal_cell_to_dual_node_map[cell_id].second;
face_id = next_face(cell_id, face_id, node_id);
cell_id = next_cell(cell_id, face_id);
} while (face_id != first_face_id);
}
});
+
+ dual_descriptor.setCellToNodeMatrix(ConnectivityMatrix(cell_to_node_row, cell_to_node_list));
}
template <>
@@ -277,39 +302,28 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
const auto& primal_node_list = primal_ref_node_list.list();
const std::vector<NodeId> dual_node_list = [&]() {
- std::vector<NodeId> dual_node_list;
+ std::vector<NodeId> tmp_dual_node_list;
for (size_t i_primal_node = 0; i_primal_node < primal_node_list.size(); ++i_primal_node) {
auto primal_node_id = primal_node_list[i_primal_node];
const auto i_dual_node =
primal_boundary_node_id_to_dual_node_id_map.find(primal_connectivity.nodeNumber()[primal_node_id]);
if (i_dual_node != primal_boundary_node_id_to_dual_node_id_map.end()) {
- dual_node_list.push_back(i_dual_node->second);
+ tmp_dual_node_list.push_back(i_dual_node->second);
}
}
- return dual_node_list;
+ return tmp_dual_node_list;
}();
if (parallel::allReduceOr(dual_node_list.size() > 0)) {
- dual_descriptor.addRefItemList(RefNodeList{primal_ref_node_list.refId(), convert_to_array(dual_node_list),
- primal_ref_node_list.isBoundary()});
+ auto dual_node_array = convert_to_array(dual_node_list);
+ dual_descriptor.addRefItemList(
+ RefNodeList{primal_ref_node_list.refId(), dual_node_array, primal_ref_node_list.isBoundary()});
}
}
}
- using Face = ConnectivityFace<2>;
-
- 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 < dual_descriptor.face_to_node_vector.size(); ++l) {
- const auto& node_vector = dual_descriptor.face_to_node_vector[l];
-
- face_to_id_map[Face(node_vector, dual_descriptor.node_number_vector)] = l;
- }
- return face_to_id_map;
- }();
-
for (size_t i_face_list = 0; i_face_list < primal_connectivity.template numberOfRefItemList<ItemType::face>();
++i_face_list) {
const auto& primal_ref_face_list = primal_connectivity.template refItemList<ItemType::face>(i_face_list);
@@ -321,21 +335,23 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
bounday_face_dual_node_id_list[i_face] = m_primal_face_to_dual_node_map[primal_face_list[i_face]].second;
}
- std::vector<bool> is_dual_node_from_boundary_face(dual_descriptor.node_number_vector.size(), false);
+ std::vector<bool> is_dual_node_from_boundary_face(dual_descriptor.nodeNumberVector().size(), false);
for (size_t i_face = 0; i_face < bounday_face_dual_node_id_list.size(); ++i_face) {
is_dual_node_from_boundary_face[bounday_face_dual_node_id_list[i_face]] = true;
}
- std::vector<bool> is_dual_node_from_boundary_node(dual_descriptor.node_number_vector.size(), false);
+ std::vector<bool> is_dual_node_from_boundary_node(dual_descriptor.nodeNumberVector().size(), false);
for (size_t i_node = 0; i_node < m_primal_boundary_node_to_dual_node_map.size(); ++i_node) {
is_dual_node_from_boundary_node[m_primal_boundary_node_to_dual_node_map[i_node].second] = true;
}
+ const auto& dual_face_to_node_matrix = dual_descriptor.faceToNodeMatrix();
+
std::vector<FaceId> dual_face_list;
dual_face_list.reserve(2 * primal_face_list.size());
- for (size_t i_dual_face = 0; i_dual_face < dual_descriptor.face_to_node_vector.size(); ++i_dual_face) {
- const NodeId dual_node_0 = dual_descriptor.face_to_node_vector[i_dual_face][0];
- const NodeId dual_node_1 = dual_descriptor.face_to_node_vector[i_dual_face][1];
+ for (size_t i_dual_face = 0; i_dual_face < dual_face_to_node_matrix.numberOfRows(); ++i_dual_face) {
+ const NodeId dual_node_0 = dual_face_to_node_matrix[i_dual_face][0];
+ const NodeId dual_node_1 = dual_face_to_node_matrix[i_dual_face][1];
if ((is_dual_node_from_boundary_face[dual_node_0] and is_dual_node_from_boundary_node[dual_node_1]) or
(is_dual_node_from_boundary_node[dual_node_0] and is_dual_node_from_boundary_face[dual_node_1])) {
@@ -352,42 +368,61 @@ MedianDualConnectivityBuilder::_buildConnectivityFrom<2>(const IConnectivity& i_
}
}
- const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
- const size_t primal_number_of_cells = primal_connectivity.numberOfCells();
+ const auto& primal_node_owner = primal_connectivity.nodeOwner();
- dual_descriptor.node_owner_vector.resize(dual_descriptor.node_number_vector.size());
+ dual_descriptor.setNodeOwnerVector([&] {
+ Array<int> node_owner_vector(dual_descriptor.nodeNumberVector().size());
- const auto& primal_node_owner = primal_connectivity.nodeOwner();
- for (NodeId primal_node_id = 0; primal_node_id < primal_connectivity.numberOfNodes(); ++primal_node_id) {
- dual_descriptor.node_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
- }
- const auto& primal_cell_owner = primal_connectivity.cellOwner();
- for (CellId primal_cell_id = 0; primal_cell_id < primal_number_of_cells; ++primal_cell_id) {
- dual_descriptor.node_owner_vector[primal_number_of_nodes + primal_cell_id] = primal_cell_owner[primal_cell_id];
- }
+ node_owner_vector.fill(-1);
+ for (size_t i_primal_node_to_dual_node = 0;
+ i_primal_node_to_dual_node < m_primal_boundary_node_to_dual_node_map.size(); ++i_primal_node_to_dual_node) {
+ const auto& [primal_node_id, dual_node_id] = m_primal_boundary_node_to_dual_node_map[i_primal_node_to_dual_node];
+ node_owner_vector[dual_node_id] = primal_node_owner[primal_node_id];
+ }
- dual_descriptor.cell_owner_vector.resize(dual_descriptor.cell_number_vector.size());
- for (NodeId primal_node_id = 0; primal_node_id < primal_number_of_nodes; ++primal_node_id) {
- dual_descriptor.cell_owner_vector[primal_node_id] = primal_node_owner[primal_node_id];
- }
+ const auto& primal_face_owner = primal_connectivity.faceOwner();
+ for (size_t i_primal_face_to_dual_node = 0; i_primal_face_to_dual_node < m_primal_face_to_dual_node_map.size();
+ ++i_primal_face_to_dual_node) {
+ const auto& [primal_face_id, dual_node_id] = m_primal_face_to_dual_node_map[i_primal_face_to_dual_node];
+ node_owner_vector[dual_node_id] = primal_face_owner[primal_face_id];
+ }
- {
- std::vector<int> face_cell_owner(dual_descriptor.face_number_vector.size());
- std::fill(std::begin(face_cell_owner), std::end(face_cell_owner), parallel::size());
+ const auto& primal_cell_owner = primal_connectivity.cellOwner();
+ for (size_t i_primal_cell_to_dual_node = 0; i_primal_cell_to_dual_node < m_primal_cell_to_dual_node_map.size();
+ ++i_primal_cell_to_dual_node) {
+ const auto& [primal_cell_id, dual_node_id] = m_primal_cell_to_dual_node_map[i_primal_cell_to_dual_node];
+ node_owner_vector[dual_node_id] = primal_cell_owner[primal_cell_id];
+ }
- for (size_t i_cell = 0; i_cell < dual_descriptor.cell_to_face_vector.size(); ++i_cell) {
- const auto& cell_face_list = dual_descriptor.cell_to_face_vector[i_cell];
- for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
- const size_t face_id = cell_face_list[i_face];
- face_cell_owner[face_id] = std::min(face_cell_owner[face_id], dual_descriptor.cell_number_vector[i_cell]);
- }
+ return node_owner_vector;
+ }());
+
+ dual_descriptor.setCellOwnerVector([&] {
+ Array<int> cell_owner_vector(dual_descriptor.cellNumberVector().size());
+ cell_owner_vector.fill(-1);
+ for (size_t i_primal_cell_to_dual_node = 0; i_primal_cell_to_dual_node < m_primal_node_to_dual_cell_map.size();
+ ++i_primal_cell_to_dual_node) {
+ const auto& [primal_node_id, dual_cell_id] = m_primal_node_to_dual_cell_map[i_primal_cell_to_dual_node];
+ cell_owner_vector[dual_cell_id] = primal_node_owner[primal_node_id];
}
+ return cell_owner_vector;
+ }());
- dual_descriptor.face_owner_vector.resize(face_cell_owner.size());
- for (size_t i_face = 0; i_face < face_cell_owner.size(); ++i_face) {
- dual_descriptor.face_owner_vector[i_face] = dual_descriptor.cell_owner_vector[face_cell_owner[i_face]];
+ dual_descriptor.setFaceOwnerVector([&] {
+ const auto& dual_cell_to_face_matrix = dual_descriptor.cellToFaceMatrix();
+ const auto& dual_cell_owner_vector = dual_descriptor.cellOwnerVector();
+
+ Array<int> face_owner_vector(dual_descriptor.faceNumberVector().size());
+ face_owner_vector.fill(parallel::size());
+ for (size_t i_cell = 0; i_cell < dual_cell_to_face_matrix.numberOfRows(); ++i_cell) {
+ const auto& cell_face_list = dual_cell_to_face_matrix[i_cell];
+ for (size_t i_face = 0; i_face < cell_face_list.size(); ++i_face) {
+ const size_t face_id = cell_face_list[i_face];
+ face_owner_vector[face_id] = std::min(face_owner_vector[face_id], dual_cell_owner_vector[i_cell]);
+ }
}
- }
+ return face_owner_vector;
+ }());
m_connectivity = ConnectivityType::build(dual_descriptor);
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index a3a8440bd25499d8aea34a08e7c836b66a4c3a32..973a6099afb49595d02403c697bbb927d2a7ba63 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -132,8 +132,14 @@ MeshBuilderBase::_checkMesh() const
if (intersection.size() > 1) {
std::ostringstream error_msg;
error_msg << "invalid mesh.\n\tFollowing faces\n";
- for (FaceId face_id : intersection) {
- error_msg << "\t - id=" << face_id << " number=" << connectivity.faceNumber()[face_id] << '\n';
+ for (FaceId intersection_face_id : intersection) {
+ error_msg << "\t - id=" << intersection_face_id
+ << " number=" << connectivity.faceNumber()[intersection_face_id] << '\n';
+ error_msg << "\t nodes:";
+ for (size_t i = 0; i < face_to_node_matrix[intersection_face_id].size(); ++i) {
+ error_msg << ' ' << face_to_node_matrix[intersection_face_id][i];
+ }
+ error_msg << '\n';
}
error_msg << "\tare duplicated";
throw NormalError(error_msg.str());
diff --git a/src/utils/Array.hpp b/src/utils/Array.hpp
index 7535b673565fa04aadba3935420f6946f5ad5544..3ee64824a650e18714a9b4b8523afb78cc89f643 100644
--- a/src/utils/Array.hpp
+++ b/src/utils/Array.hpp
@@ -24,6 +24,18 @@ class [[nodiscard]] Array
const size_t m_size;
public:
+ friend std::ostream&
+ operator<<(std::ostream& os, const UnsafeArrayView& x)
+ {
+ if (x.size() > 0) {
+ os << 0 << ':' << NaNHelper(x[0]);
+ }
+ for (size_t i = 1; i < x.size(); ++i) {
+ os << ' ' << i << ':' << NaNHelper(x[i]);
+ }
+ return os;
+ }
+
[[nodiscard]] PUGS_INLINE size_t
size() const
{
diff --git a/src/utils/Table.hpp b/src/utils/Table.hpp
index 6972ea1030dedbdd30c039e0880a3773024603dc..a4a134dfd1d5838766b062eb429c725eb66bbb33 100644
--- a/src/utils/Table.hpp
+++ b/src/utils/Table.hpp
@@ -63,6 +63,19 @@ class [[nodiscard]] Table
Assert(row < table.numberOfRows(), "required row view is not contained in the Table");
}
+ friend std::ostream&
+ operator<<(std::ostream& os, const UnsafeRowView& x)
+ {
+ if (x.size() > 0) {
+ os << 0 << ':' << NaNHelper(x[0]);
+ }
+ for (size_t i = 1; i < x.size(); ++i) {
+ os << ' ' << i << ':' << NaNHelper(x[i]);
+ }
+
+ return os;
+ }
+
// To try to keep these views close to the initial array one
// forbids copy constructor and take benefits of C++-17 copy
// elisions.
@@ -113,6 +126,19 @@ class [[nodiscard]] Table
}
}
+ friend std::ostream&
+ operator<<(std::ostream& os, const RowView& x)
+ {
+ if (x.size() > 0) {
+ os << 0 << ':' << NaNHelper(x[0]);
+ }
+ for (size_t i = 1; i < x.size(); ++i) {
+ os << ' ' << i << ':' << NaNHelper(x[i]);
+ }
+
+ return os;
+ }
+
RowView(const UnsafeTableView& table_view, index_type row) : m_table_view{table_view}, m_row{row}
{
Assert(row < m_table_view.numberOfRows(), "required row view is not contained in the Table view");
@@ -155,6 +181,19 @@ class [[nodiscard]] Table
return m_table(m_row_begin + i, m_column_begin + j);
}
+ friend std::ostream&
+ operator<<(std::ostream& os, const UnsafeTableView& t)
+ {
+ for (size_t i = 0; i < t.numberOfRows(); ++i) {
+ os << i << '|';
+ for (size_t j = 0; j < t.numberOfColumns(); ++j) {
+ os << ' ' << j << ':' << NaNHelper(t(i, j));
+ }
+ os << '\n';
+ }
+ return os;
+ }
+
PUGS_INLINE void
fill(const DataType& data) const
{
diff --git a/tests/test_BuiltinFunctionProcessor.cpp b/tests/test_BuiltinFunctionProcessor.cpp
index 0168633cea1d2d2e764b6b835093d1f17504024d..8de9ec2c2ed9e61479aa2d7165bf6e506fd6caf6 100644
--- a/tests/test_BuiltinFunctionProcessor.cpp
+++ b/tests/test_BuiltinFunctionProcessor.cpp
@@ -385,15 +385,15 @@ let d:R, d = det(A);
tested_function_set.insert("det:R^3x3");
std::string_view data = R"(
import math;
-let A:R^3x3, A = [[-2.5, 2.9,-1.3],
- [ 3.2, 2.3, 2.7],
- [-2.6, 5.2,-3.5]];
+let A:R^3x3, A = [[-2, 2,-1],
+ [ 3, 2, 2],
+ [-2, 5,-3]];
let d:R, d = det(A);
)";
CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "d",
- det(TinyMatrix<3>{-2.5, 2.9, -1.3, //
- +3.2, 2.3, +2.7, //
- -2.6, 5.2, -3.5}));
+ det(TinyMatrix<3>{-2, 2, -1, //
+ +3, 2, +2, //
+ -2, 5, -3}));
}
{ // trace
@@ -461,15 +461,15 @@ let invA:R^2x2, invA = inverse(A);
tested_function_set.insert("inverse:R^3x3");
std::string_view data = R"(
import math;
-let A:R^3x3, A = [[-2.5, 2.9,-1.3],
- [ 3.2, 2.3, 2.7],
- [-2.6, 5.2,-3.5]];
+let A:R^3x3, A = [[-2, 2,-1],
+ [ 3, 2, 2],
+ [-2, 5,-3]];
let invA:R^3x3, invA = inverse(A);
)";
CHECK_BUILTIN_FUNCTION_EVALUATION_RESULT(data, "invA",
- inverse(TinyMatrix<3>{-2.5, 2.9, -1.3, //
- +3.2, 2.3, +2.7, //
- -2.6, 5.2, -3.5}));
+ inverse(TinyMatrix<3>{-2, 2, -1, //
+ +3, 2, +2, //
+ -2, 5, -3}));
}
{ // transpose
diff --git a/tests/test_Connectivity.cpp b/tests/test_Connectivity.cpp
index 60e894a83a367b6d07ef80939ac0860095dfd322..b0465d08489c644685500d59b145e12594966522 100644
--- a/tests/test_Connectivity.cpp
+++ b/tests/test_Connectivity.cpp
@@ -3,6 +3,7 @@
#include <MeshDataBaseForTests.hpp>
#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
#include <mesh/ItemValue.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
@@ -670,6 +671,318 @@ TEST_CASE("Connectivity", "[mesh]")
}
}
+ SECTION("item ordering")
+ {
+ SECTION("1D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all1DMeshes();
+
+ for (const auto& named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ SECTION("cell -> nodes")
+ {
+ auto mesh = named_mesh.mesh();
+ auto xr = mesh->xr();
+
+ const Connectivity<1>& connectivity = mesh->connectivity();
+
+ auto cell_to_node_matrix = connectivity.cellToNodeMatrix();
+
+ bool is_correct = true;
+
+ for (CellId cell_id = 0; cell_id < connectivity.numberOfCells(); ++cell_id) {
+ if (xr[cell_to_node_matrix[cell_id][1]][0] < xr[cell_to_node_matrix[cell_id][0]][0]) {
+ is_correct = false;
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<1>& connectivity = mesh->connectivity();
+
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto cell_node_list = node_to_cell_matrix[node_id];
+ for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+ is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+ }
+ }
+ }
+
+ SECTION("2D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all2DMeshes();
+
+ for (const auto& named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ SECTION("face -> nodes")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh->connectivity();
+
+ auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+ auto node_number = connectivity.nodeNumber();
+
+ bool is_correct = true;
+
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ if (node_number[face_node_list[1]] < node_number[face_node_list[0]]) {
+ is_correct = false;
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> faces")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh->connectivity();
+
+ auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+ auto face_number = connectivity.faceNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto face_node_list = node_to_face_matrix[node_id];
+ for (size_t i_node = 0; i_node < face_node_list.size() - 1; ++i_node) {
+ is_correct &= (face_number[face_node_list[i_node]] < face_number[face_node_list[i_node + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh->connectivity();
+
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto cell_node_list = node_to_cell_matrix[node_id];
+ for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+ is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("face -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<2>& connectivity = mesh->connectivity();
+
+ auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto cell_face_list = face_to_cell_matrix[face_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size() - 1; ++i_face) {
+ is_correct &= (cell_number[cell_face_list[i_face]] < cell_number[cell_face_list[i_face + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+ }
+ }
+ }
+
+ SECTION("3D")
+ {
+ std::array mesh_list = MeshDataBaseForTests::get().all3DMeshes();
+
+ for (const auto& named_mesh : mesh_list) {
+ SECTION(named_mesh.name())
+ {
+ SECTION("edge -> nodes")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto edge_to_node_matrix = connectivity.edgeToNodeMatrix();
+ auto node_number = connectivity.nodeNumber();
+
+ bool is_correct = true;
+
+ for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+ auto edge_node_list = edge_to_node_matrix[edge_id];
+ if (node_number[edge_node_list[1]] < node_number[edge_node_list[0]]) {
+ is_correct = false;
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("face -> nodes")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto face_to_node_matrix = connectivity.faceToNodeMatrix();
+ auto node_number = connectivity.nodeNumber();
+
+ bool is_correct = true;
+
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto face_node_list = face_to_node_matrix[face_id];
+ for (size_t i = 1; i < face_node_list.size() - 1; ++i) {
+ if (node_number[face_node_list[i]] < node_number[face_node_list[0]]) {
+ is_correct = false;
+ }
+ }
+ for (size_t i = 2; i < face_node_list.size() - 1; ++i) {
+ if (node_number[face_node_list[i]] < node_number[face_node_list[1]]) {
+ is_correct = false;
+ }
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> edges")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ REQUIRE(checkConnectivityOrdering(connectivity));
+
+ auto node_to_edge_matrix = connectivity.nodeToEdgeMatrix();
+ auto edge_number = connectivity.edgeNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto edge_node_list = node_to_edge_matrix[node_id];
+ for (size_t i_node = 0; i_node < edge_node_list.size() - 1; ++i_node) {
+ is_correct &= (edge_number[edge_node_list[i_node]] < edge_number[edge_node_list[i_node + 1]]);
+ }
+ }
+
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> faces")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto node_to_face_matrix = connectivity.nodeToFaceMatrix();
+ auto face_number = connectivity.faceNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto face_node_list = node_to_face_matrix[node_id];
+ for (size_t i_node = 0; i_node < face_node_list.size() - 1; ++i_node) {
+ is_correct &= (face_number[face_node_list[i_node]] < face_number[face_node_list[i_node + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("node -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto node_to_cell_matrix = connectivity.nodeToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (NodeId node_id = 0; node_id < connectivity.numberOfNodes(); ++node_id) {
+ auto cell_node_list = node_to_cell_matrix[node_id];
+ for (size_t i_node = 0; i_node < cell_node_list.size() - 1; ++i_node) {
+ is_correct &= (cell_number[cell_node_list[i_node]] < cell_number[cell_node_list[i_node + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("edge -> faces")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto edge_to_face_matrix = connectivity.edgeToFaceMatrix();
+ auto face_number = connectivity.faceNumber();
+
+ bool is_correct = true;
+ for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+ auto face_edge_list = edge_to_face_matrix[edge_id];
+ for (size_t i_edge = 0; i_edge < face_edge_list.size() - 1; ++i_edge) {
+ is_correct &= (face_number[face_edge_list[i_edge]] < face_number[face_edge_list[i_edge + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("edge -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto edge_to_cell_matrix = connectivity.edgeToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (EdgeId edge_id = 0; edge_id < connectivity.numberOfEdges(); ++edge_id) {
+ auto cell_edge_list = edge_to_cell_matrix[edge_id];
+ for (size_t i_edge = 0; i_edge < cell_edge_list.size() - 1; ++i_edge) {
+ is_correct &= (cell_number[cell_edge_list[i_edge]] < cell_number[cell_edge_list[i_edge + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+
+ SECTION("face -> cells")
+ {
+ auto mesh = named_mesh.mesh();
+
+ const Connectivity<3>& connectivity = mesh->connectivity();
+
+ auto face_to_cell_matrix = connectivity.faceToCellMatrix();
+ auto cell_number = connectivity.cellNumber();
+
+ bool is_correct = true;
+ for (FaceId face_id = 0; face_id < connectivity.numberOfFaces(); ++face_id) {
+ auto cell_face_list = face_to_cell_matrix[face_id];
+ for (size_t i_face = 0; i_face < cell_face_list.size() - 1; ++i_face) {
+ is_correct &= (cell_number[cell_face_list[i_face]] < cell_number[cell_face_list[i_face + 1]]);
+ }
+ }
+ REQUIRE(is_correct);
+ }
+ }
+ }
+ }
+ }
+
SECTION("ItemLocalNumbersInTheirSubItems")
{
auto check_item_local_numbers_in_their_subitems = [](auto item_to_subitem_matrix, auto subitem_to_item_matrix,
diff --git a/tests/test_DiamondDualConnectivityBuilder.cpp b/tests/test_DiamondDualConnectivityBuilder.cpp
index 70b8f75440061530f386654a4bbf9e3f5342d33e..45d9cd700ca7b5ae73e6c6fd22db024e7de3126c 100644
--- a/tests/test_DiamondDualConnectivityBuilder.cpp
+++ b/tests/test_DiamondDualConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
@@ -54,6 +55,8 @@ TEST_CASE("DiamondDualConnectivityBuilder", "[mesh]")
REQUIRE(dual_connectivity.numberOfFaces() == 220);
REQUIRE(dual_connectivity.numberOfCells() == 110);
+ REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
SECTION("ref node list")
{
REQUIRE(primal_connectivity.numberOfRefItemList<ItemType::node>() == 4);
@@ -159,6 +162,8 @@ TEST_CASE("DiamondDualConnectivityBuilder", "[mesh]")
DualConnectivityManager::instance().getDiamondDualConnectivity(primal_connectivity);
const ConnectivityType& dual_connectivity = *p_diamond_dual_connectivity;
+ REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
REQUIRE(dual_connectivity.numberOfNodes() == 331);
REQUIRE(dual_connectivity.numberOfEdges() == 1461);
REQUIRE(dual_connectivity.numberOfFaces() == 1651);
diff --git a/tests/test_Dual1DConnectivityBuilder.cpp b/tests/test_Dual1DConnectivityBuilder.cpp
index 2b194a92baa821d6da4bac1e6277aaea6c1ddb1b..b7a35abfc577d93f729a0befc166c16c04f47fd3 100644
--- a/tests/test_Dual1DConnectivityBuilder.cpp
+++ b/tests/test_Dual1DConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
@@ -48,6 +49,8 @@ TEST_CASE("Dual1DConnectivityBuilder", "[mesh]")
DualConnectivityManager::instance().getDual1DConnectivity(primal_connectivity);
const ConnectivityType& dual_connectivity = *p_dual_1d_connectivity;
+ REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
REQUIRE(dual_connectivity.numberOfNodes() == 36);
REQUIRE(dual_connectivity.numberOfCells() == 35);
diff --git a/tests/test_MedianDualConnectivityBuilder.cpp b/tests/test_MedianDualConnectivityBuilder.cpp
index a95e6ad22f8a9b2605dec5c28e7c7704d40aed6f..6d091bd037e4ac57635757b98671eff9fe6eb643 100644
--- a/tests/test_MedianDualConnectivityBuilder.cpp
+++ b/tests/test_MedianDualConnectivityBuilder.cpp
@@ -5,6 +5,7 @@
#include <mesh/DualConnectivityManager.hpp>
#include <mesh/Connectivity.hpp>
+#include <mesh/ConnectivityUtils.hpp>
#include <mesh/ItemValueUtils.hpp>
#include <mesh/Mesh.hpp>
@@ -50,6 +51,8 @@ TEST_CASE("MedianDualConnectivityBuilder", "[mesh]")
DualConnectivityManager::instance().getMedianDualConnectivity(primal_connectivity);
const ConnectivityType& dual_connectivity = *p_median_dual_connectivity;
+ REQUIRE(checkConnectivityOrdering(dual_connectivity));
+
REQUIRE(dual_connectivity.numberOfNodes() == 192);
REQUIRE(dual_connectivity.numberOfFaces() == 244);
REQUIRE(dual_connectivity.numberOfCells() == 53);