diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index 9f6ea50cab88390669e3471ad3c70f25e870e390..79fb0a2011bceae67a20e7b877f2896a2af4e66e 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -5,6 +5,8 @@
#include <utils/Array.hpp>
#include <utils/Messenger.hpp>
+#include <utils/Timer.hpp>
+
template <size_t Dimension>
NodeValue<TinyVector<Dimension>>
CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<Dimension>&,
@@ -67,6 +69,9 @@ CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
const TinyVector<3, uint64_t>& cell_size,
const IConnectivity& connectivity) const
{
+ Timer t;
+ t.start();
+
const TinyVector<3, uint64_t> node_size = [&] {
TinyVector node_size{cell_size};
for (size_t i = 0; i < 3; ++i) {
@@ -95,6 +100,8 @@ CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
}
});
+ std::cout << "CartesianMeshBuilder::_getNodeCoordinates t = " << t << '\n';
+
return xr;
}
@@ -104,6 +111,8 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
const TinyVector<Dimension>& b,
const TinyVector<Dimension, uint64_t>& cell_size)
{
+ Timer t;
+ t.start();
static_assert(Dimension <= 3, "unexpected dimension");
LogicalConnectivityBuilder logical_connectivity_builder{cell_size};
@@ -117,6 +126,8 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
NodeValue<TinyVector<Dimension>> xr = _getNodeCoordinates(a, b, cell_size, connectivity);
m_mesh = std::make_shared<Mesh<ConnectivityType>>(p_connectivity, xr);
+
+ std::cout << "_buildCartesianMesh t = " << t << '\n';
}
template <size_t Dimension>
diff --git a/src/mesh/ConnectivityBuilderBase.cpp b/src/mesh/ConnectivityBuilderBase.cpp
index 55c56b3a4b78f28ecdf8613d04d4d0348b9e5c4a..9a02e30931998e629b1f6b480ff28191d6453405 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,538 @@ 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;
+
+ size_t total_number_of_faces = 0;
+
for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
const auto& cell_nodes = descriptor.cell_to_node_vector[j];
if constexpr (Dimension == 2) {
+ total_number_of_faces += cell_nodes.size();
+ } else if constexpr (Dimension == 3) {
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()));
+ 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(descriptor.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()));
+ }
+ }
+
+ 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 < descriptor.cell_to_node_vector.size(); ++j) {
+ switch (descriptor.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 = descriptor.cell_to_node_vector[j];
+ count_number_of_face_by_node += 1 * cell_nodes.size() + cell_nodes.size() * 3;
+ break;
+ }
+ case CellType::Diamond: {
+ const auto& cell_nodes = descriptor.cell_to_node_vector[j];
+ count_number_of_face_by_node += cell_nodes.size() * 3 * 2;
+ 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());
+ }
+ }
+ }
+ return count_number_of_face_by_node;
+ }
+ }();
- // 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()));
+ Array<unsigned int> faces_node_list(total_number_of_face_by_node);
+ Array<unsigned int> face_ending(total_number_of_faces + 1);
+ size_t i_face = 0;
+ face_ending[0] = 0;
- // 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 short> cell_nb_faces(descriptor.cell_to_node_vector.size());
+ {
+ size_t i_face_node = 0;
+ for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
+ const auto& cell_nodes = descriptor.cell_to_node_vector[j];
- // 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()));
+ if constexpr (Dimension == 2) {
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Triangle: {
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ cell_nb_faces[j] = 3;
+ break;
+ }
+ case CellType::Quadrangle: {
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ face_ending[i_face + 1] = face_ending[i_face] + 2;
+ i_face++;
+ cell_nb_faces[j] = 4;
+ break;
+ }
+ case CellType::Polygon: {
+ for (size_t i = 0; i < cell_nodes.size(); ++i) {
+ faces_node_list[i_face_node] = cell_nodes[i];
+ i_face_node++;
+ faces_node_list[i_face_node] = cell_nodes[(i + 1) % cell_nodes.size()];
+ i_face_node++;
+
+ face_ending[i_face + 1] = face_ending[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(descriptor.cell_type_vector[j]) << ": unexpected cell type in dimension 2";
+ throw UnexpectedError(error_msg.str());
+ }
+ }
+ } else if constexpr (Dimension == 3) {
+ switch (descriptor.cell_type_vector[j]) {
+ case CellType::Hexahedron: {
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[4];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+ faces_node_list[i_face_node++] = cell_nodes[6];
+ faces_node_list[i_face_node++] = cell_nodes[7];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[4];
+ faces_node_list[i_face_node++] = cell_nodes[7];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[6];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+ faces_node_list[i_face_node++] = cell_nodes[4];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[7];
+ faces_node_list[i_face_node++] = cell_nodes[6];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ cell_nb_faces[j] = 6;
+ break;
+ }
+ case CellType::Tetrahedron: {
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_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()));
+ face_ending[i_face + 1] = face_ending[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()));
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_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()));
+ face_ending[i_face + 1] = face_ending[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()));
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_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()));
+ face_ending[i_face + 1] = face_ending[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()));
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_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];
- }
+ face_ending[i_face + 1] = face_ending[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: {
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 3;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[4];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 3;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+ faces_node_list[i_face_node++] = cell_nodes[4];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[1];
+ faces_node_list[i_face_node++] = cell_nodes[4];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+
+ face_ending[i_face + 1] = face_ending[i_face] + 4;
+ i_face++;
+
+ faces_node_list[i_face_node++] = cell_nodes[2];
+ faces_node_list[i_face_node++] = cell_nodes[0];
+ faces_node_list[i_face_node++] = cell_nodes[3];
+ faces_node_list[i_face_node++] = cell_nodes[5];
+
+ face_ending[i_face + 1] = face_ending[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) {
+ faces_node_list[i_face_node++] = base_nodes[i];
+ }
+
+ face_ending[i_face + 1] = face_ending[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()));
+ faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
+ faces_node_list[i_face_node++] = base_nodes[i_node];
+ faces_node_list[i_face_node++] = pyramid_vertex;
+
+ face_ending[i_face + 1] = face_ending[i_face] + 3;
+ i_face++;
}
+ break;
}
+ case CellType::Diamond: {
+ std::vector<unsigned int> base_nodes;
+ std::copy_n(cell_nodes.begin() + 1, cell_nodes.size() - 2, std::back_inserter(base_nodes));
+
+ { // top faces
+ const auto top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_node = 0; i_node < base_nodes.size(); ++i_node) {
+ faces_node_list[i_face_node++] = base_nodes[i_node];
+ faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
+ faces_node_list[i_face_node++] = top_vertex;
+
+ face_ending[i_face + 1] = face_ending[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 < base_nodes.size(); ++i_node) {
+ faces_node_list[i_face_node++] = base_nodes[(i_node + 1) % base_nodes.size()];
+ faces_node_list[i_face_node++] = base_nodes[i_node];
+ faces_node_list[i_face_node++] = bottom_vertex;
+
+ face_ending[i_face + 1] = face_ending[i_face] + 3;
+ i_face++;
+ }
}
+ cell_nb_faces[j] = 2 * (cell_nodes.size() - 2);
+ 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());
+ }
}
- 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> face_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[faces_node_list[2 * i_face]] > node_number_vector[faces_node_list[2 * i_face + 1]]) {
+ std::swap(faces_node_list[2 * i_face], faces_node_list[2 * i_face + 1]);
+ face_reversed[i_face] = true;
+ } else {
+ face_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 = face_ending[i_face + 1] - face_ending[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[faces_node_list[face_ending[i_face] + i_face_node]] <
+ node_number_vector[faces_node_list[face_ending[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] = faces_node_list[face_ending[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] =
+ faces_node_list[face_ending[i_face] + i_node];
+ }
+
+ for (size_t i_node = 0; i_node < face_node_number; ++i_node) {
+ faces_node_list[face_ending[i_face] + i_node] = buffer[i_node];
+ }
+ }
+
+ if (node_number_vector[faces_node_list[face_ending[i_face] + 1]] >
+ node_number_vector[faces_node_list[face_ending[i_face + 1] - 1]]) {
+ for (size_t i_node = 1; i_node <= (face_node_number + 1) / 2 - 1; ++i_node) {
+ std::swap(faces_node_list[face_ending[i_face] + i_node], faces_node_list[face_ending[i_face + 1] - i_node]);
+ }
+
+ face_reversed[i_face] = true;
+ } else {
+ face_reversed[i_face] = false;
}
}
}
- {
- 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]);
- }
- 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> 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 = face_ending[i_face]; i_node_face < face_ending[i_face + 1]; ++i_node_face) {
+ node_to_duplicate_face_row[faces_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_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]);
+ Array<unsigned int> node_face_row_idx(node_number_vector.size());
+ node_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 = face_ending[i_face]; i_node_face < face_ending[i_face + 1]; ++i_node_face) {
+ const size_t node_id = faces_node_list[i_node_face];
+ node_duplicate_face_list[node_to_duplicate_face_row[node_id] + node_face_row_idx[node_id]] = i_face;
+ node_face_row_idx[node_id]++;
+ }
}
- 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> 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 ((face_ending[i_face + 1] - face_ending[i_face]) != (face_ending[j_face + 1] - face_ending[j_face])) {
+ return false;
+ } else {
+ auto i_face_node = face_ending[i_face];
+ auto j_face_node = face_ending[j_face];
+ while (i_face_node < face_ending[i_face + 1]) {
+ if (faces_node_list[i_face_node] != faces_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.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;
+ 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];
+ }
+ }
+ }
+
+ Array<unsigned int> cell_face_begin(cell_nb_faces.size() + 1);
+ cell_face_begin[0] = 0;
+ for (size_t i = 0; i < cell_nb_faces.size(); ++i) {
+ cell_face_begin[i + 1] = cell_face_begin[i] + cell_nb_faces[i];
+ }
+
+ 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]);
+
+ for (size_t i_face = 0; i_face < cell_nb_faces[j]; ++i_face) {
+ descriptor.cell_to_face_vector[j][i_face] = dup_face_to_face_id[cell_face_begin[j] + i_face];
}
}
{
// 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.resize(total_number_of_faces - nb_dup_faces);
+ for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
descriptor.face_number_vector[l] = l;
}
}
+
+ {
+ descriptor.face_to_node_vector.resize(total_number_of_faces - nb_dup_faces);
+ int l = 0;
+
+ for (size_t i_face = 0; i_face < total_number_of_faces; ++i_face) {
+ if (dup_face_to_face[i_face] == i_face) {
+ descriptor.face_to_node_vector[l].resize(face_ending[i_face + 1] - face_ending[i_face]);
+
+ for (size_t i_face_node = 0; i_face_node < face_ending[i_face + 1] - face_ending[i_face]; ++i_face_node) {
+ descriptor.face_to_node_vector[l][i_face_node] = faces_node_list[face_ending[i_face] + i_face_node];
+ }
+ l++;
+ }
+ }
+ }
+
+ {
+ descriptor.cell_face_is_reversed_vector.resize(descriptor.cell_to_node_vector.size());
+ size_t l = 0;
+ for (CellId j = 0; j < descriptor.cell_face_is_reversed_vector.size(); ++j) {
+ descriptor.cell_face_is_reversed_vector[j].resize(cell_nb_faces[j]);
+ for (size_t i_face = 0; i_face < cell_nb_faces[j]; ++i_face) {
+ descriptor.cell_face_is_reversed_vector[j][i_face] = face_reversed[l];
+ ++l;
+ }
+ }
+ }
}
template <size_t Dimension>
@@ -257,76 +553,248 @@ 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()));
+ {
+ 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();
}
}
- std::unordered_map<Edge, EdgeId, typename Edge::Hash> edge_id_map;
+ Array<unsigned int> face_edge_ending(face_nb_edges.size() + 1);
+ face_edge_ending[0] = 0;
+ for (size_t i_face = 0; i_face < face_nb_edges.size(); ++i_face) {
+ face_edge_ending[i_face + 1] = face_edge_ending[i_face] + face_nb_edges[i_face];
+ }
+
+ const size_t total_number_of_face_edges = face_edge_ending[face_edge_ending.size() - 1];
+
+ const size_t total_number_of_node_by_face_edges = 2 * total_number_of_face_edges;
+
+ Array<unsigned int> face_edges_node_list(total_number_of_node_by_face_edges);
{
- 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]);
+ size_t i_edge_node = 0;
+ for (size_t i_face = 0; i_face < face_nb_edges.size(); ++i_face) {
+ for (size_t i_node = 0; i_node < descriptor.face_to_node_vector[i_face].size() - 1; ++i_node) {
+ face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][i_node];
+ face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][i_node + 1];
+ }
+ face_edges_node_list[i_edge_node++] =
+ descriptor.face_to_node_vector[i_face][descriptor.face_to_node_vector[i_face].size() - 1];
+ face_edges_node_list[i_edge_node++] = descriptor.face_to_node_vector[i_face][0];
}
- 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;
+ }
+
+ 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 (descriptor.node_number_vector[face_edges_node_list[2 * i_edge]] >
+ descriptor.node_number_vector[face_edges_node_list[2 * i_edge + 1]]) {
+ std::swap(face_edges_node_list[2 * i_edge], face_edges_node_list[2 * i_edge + 1]);
+ face_edge_is_reversed[i_edge] = true;
+ } else {
+ face_edge_is_reversed[i_edge] = false;
}
}
+ const size_t total_number_of_edges = face_edge_is_reversed.size();
+ Array<unsigned int> node_to_face_edge_row(descriptor.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_edges_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(descriptor.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_edges_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_edges_node_list[i_edge_first_node] == face_edges_node_list[j_edge_first_node]) and
+ (face_edges_node_list[i_edge_first_node + 1] == face_edges_node_list[j_edge_first_node + 1]));
+ };
+
+ const auto& node_number_vector = descriptor.node_number_vector;
+ 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) {
+ for (size_t j_node_face = i_node_face + 1; j_node_face < node_to_face_edge_row[i_node + 1]; ++j_node_face) {
+ if (face_edge_to_edge[node_to_face_edge_list[i_node_face]] ==
+ face_edge_to_edge[node_to_face_edge_list[j_node_face]])
+ continue;
+ if (is_same_face(node_to_face_edge_list[i_node_face], node_to_face_edge_list[j_node_face])) {
+ face_edge_to_edge[node_to_face_edge_list[j_node_face]] =
+ face_edge_to_edge[node_to_face_edge_list[i_node_face]];
+ 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<NodeId> edge_to_node_list(2 * (total_number_of_edges - nb_duplicate_edges));
+ {
+ 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_edges_node_list[2 * i_dup_edge];
+ edge_to_node_list[2 * edge_id + 1] = face_edges_node_list[2 * i_dup_edge + 1];
+ }
}
}
{
// Edge numbers may change if numbers are provided in the file
- descriptor.edge_number_vector.resize(edge_faces_map.size());
- for (size_t e = 0; e < edge_faces_map.size(); ++e) {
- descriptor.edge_number_vector[e] = e;
+ descriptor.edge_number_vector.resize(total_number_of_edges - nb_duplicate_edges);
+ for (size_t i_edge = 0; i_edge < descriptor.edge_number_vector.size(); ++i_edge) {
+ descriptor.edge_number_vector[i_edge] = i_edge;
+ }
+ }
+
+ descriptor.edge_to_node_vector.resize(total_number_of_edges - nb_duplicate_edges);
+ {
+ for (EdgeId edge_id = 0; edge_id < total_number_of_edges - nb_duplicate_edges; ++edge_id) {
+ descriptor.edge_to_node_vector[edge_id].resize(2);
+
+ for (size_t i_edge_node = 0; i_edge_node < 2; ++i_edge_node) {
+ descriptor.edge_to_node_vector[edge_id][i_edge_node] = edge_to_node_list[2 * edge_id + i_edge_node];
+ }
+ }
+ }
+
+ // 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.face_to_edge_vector.resize(descriptor.face_to_node_vector.size());
+ {
+ size_t i_edge = 0;
+ for (FaceId face_id = 0; face_id < descriptor.face_to_edge_vector.size(); ++face_id) {
+ descriptor.face_to_edge_vector[face_id].resize(descriptor.face_to_node_vector[face_id].size());
+
+ for (size_t i_face_edge = 0; i_face_edge < descriptor.face_to_edge_vector[face_id].size(); ++i_face_edge) {
+ descriptor.face_to_edge_vector[face_id][i_face_edge] = face_edge_to_edge[i_edge++];
+ }
}
}
+ descriptor.face_edge_is_reversed_vector.resize(descriptor.face_to_node_vector.size());
+ {
+ size_t i_edge = 0;
+ for (FaceId face_id = 0; face_id < descriptor.face_edge_is_reversed_vector.size(); ++face_id) {
+ descriptor.face_edge_is_reversed_vector[face_id].resize(descriptor.face_to_node_vector[face_id].size());
+
+ for (size_t i_face_edge = 0; i_face_edge < descriptor.face_edge_is_reversed_vector[face_id].size();
+ ++i_face_edge) {
+ descriptor.face_edge_is_reversed_vector[face_id][i_face_edge] = face_edge_is_reversed[i_edge++];
+ }
+ }
+ }
+
+ 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(descriptor.node_number_vector.size());
+ node_nb_edges.fill(0);
+
+ for (size_t node_id = 0; node_id < descriptor.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++;
+ }
+ }
+ }
+ }
+ 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(descriptor.node_number_vector.size() + 1);
+ node_to_edge_row[0] = 0;
+ for (size_t node_id = 0; node_id < descriptor.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<EdgeId> 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];
+ }
+ }
+ }
+
+ auto find_edge = [=](const NodeId& node0_id, const NodeId& node1_id) -> EdgeId {
+ auto [first_node_id, second_node_id] = [&] {
+ if (descriptor.node_number_vector[node0_id] < descriptor.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");
+ };
+
{
descriptor.cell_to_edge_vector.reserve(descriptor.cell_to_node_vector.size());
for (CellId j = 0; j < descriptor.cell_to_node_vector.size(); ++j) {
@@ -338,13 +806,9 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
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");
- }
- cell_edge_vector.push_back(i->second);
+ const auto& e = local_edge[i_edge];
+ const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+ cell_edge_vector.push_back(edge_id);
}
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
break;
@@ -355,29 +819,21 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
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);
+ const auto& e = local_edge[i_edge];
+ const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+ cell_edge_vector.push_back(edge_id);
}
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
break;
}
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}};
+ constexpr int local_edge[9][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);
+ const auto& e = local_edge[i_edge];
+ const EdgeId edge_id = find_edge(cell_nodes[e[0]], cell_nodes[e[1]]);
+ cell_edge_vector.push_back(edge_id);
}
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
break;
@@ -389,23 +845,16 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
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");
- }
- cell_edge_vector.push_back(i->second);
+ const EdgeId edge_id = find_edge(base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]);
+ cell_edge_vector.push_back(edge_id);
}
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");
- }
- cell_edge_vector.push_back(i->second);
+ const EdgeId edge_id = find_edge(base_nodes[i_edge], top_vertex);
+ cell_edge_vector.push_back(edge_id);
}
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
break;
@@ -417,37 +866,29 @@ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Co
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");
- }
- cell_edge_vector.push_back(i->second);
+ const EdgeId edge_id = find_edge(base_nodes[i_edge], base_nodes[(i_edge + 1) % base_nodes.size()]);
+ cell_edge_vector.push_back(edge_id);
}
- 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");
+ {
+ const unsigned int top_vertex = cell_nodes[cell_nodes.size() - 1];
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ const EdgeId edge_id = find_edge(base_nodes[i_edge], top_vertex);
+ cell_edge_vector.push_back(edge_id);
}
- 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 bottom_vertex = cell_nodes[0];
+ for (size_t i_edge = 0; i_edge < base_nodes.size(); ++i_edge) {
+ const EdgeId edge_id = find_edge(base_nodes[i_edge], bottom_vertex);
+ cell_edge_vector.push_back(edge_id);
}
- cell_edge_vector.push_back(i->second);
}
descriptor.cell_to_edge_vector.emplace_back(cell_edge_vector);
-
break;
}
default: {
diff --git a/src/mesh/ConnectivityDescriptor.hpp b/src/mesh/ConnectivityDescriptor.hpp
index 8f26988eadaa01dfec8aa0f14404a82819273475..0c7f60b38b9b83e95a244fb5c7d4c4b1d13bfeaf 100644
--- a/src/mesh/ConnectivityDescriptor.hpp
+++ b/src/mesh/ConnectivityDescriptor.hpp
@@ -47,8 +47,8 @@ class ConnectivityDescriptor
}
}
- std::vector<Array<bool>> cell_face_is_reversed_vector;
- std::vector<Array<bool>> face_edge_is_reversed_vector;
+ std::vector<std::vector<bool>> cell_face_is_reversed_vector;
+ std::vector<std::vector<bool>> face_edge_is_reversed_vector;
std::vector<CellType> cell_type_vector;
diff --git a/src/mesh/ConnectivityDispatcher.cpp b/src/mesh/ConnectivityDispatcher.cpp
index fd5b58c6daaf5fd44c0cb440e7af472b9f9aebbf..2cd8965f2929a8453d4f53b5d692c981fcf44906 100644
--- a/src/mesh/ConnectivityDispatcher.cpp
+++ b/src/mesh/ConnectivityDispatcher.cpp
@@ -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)
+ std::vector<std::vector<std::remove_const_t<DataType>>>& gathered_vector)
{
using MutableDataType = std::remove_const_t<DataType>;
@@ -206,7 +206,7 @@ ConnectivityDispatcher<Dimension>::_gatherFrom(
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]);
+ std::vector<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++];
}
diff --git a/src/mesh/ConnectivityDispatcher.hpp b/src/mesh/ConnectivityDispatcher.hpp
index 73ebc548d98a01e4999611927272206df4af1c23..ef1f385bcd4676a0c1afe1c9a6a53989830712bf 100644
--- a/src/mesh/ConnectivityDispatcher.hpp
+++ b/src/mesh/ConnectivityDispatcher.hpp
@@ -184,7 +184,7 @@ class ConnectivityDispatcher
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);
+ std::vector<std::vector<std::remove_const_t<DataType>>>& gathered_vector);
template <typename SubItemOfItemT>
void _buildNumberOfSubItemPerItemToRecvByProc();
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index e555ca224c1d02f11e11696542774f2e470b6432..c033433f807b2313c2dab4132f9ecc287ad413fc 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -561,8 +561,8 @@ GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData&
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;
+ for (size_t edge_id : descriptor.face_to_edge_vector[i_face]) {
+ is_boundary_edge[edge_id] = true;
}
}
}
diff --git a/src/mesh/LogicalConnectivityBuilder.cpp b/src/mesh/LogicalConnectivityBuilder.cpp
index 41092faaf90a9a3f9a24ac1c561304373ba748d0..5314f956c27de580ee227f8f8519648d4b105629 100644
--- a/src/mesh/LogicalConnectivityBuilder.cpp
+++ b/src/mesh/LogicalConnectivityBuilder.cpp
@@ -6,6 +6,8 @@
#include <utils/Array.hpp>
#include <utils/Messenger.hpp>
+#include <utils/Timer.hpp>
+
#include <unordered_map>
template <size_t Dimension>
@@ -535,6 +537,9 @@ template <>
void
LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
{
+ Timer t;
+ t.start();
+
constexpr size_t Dimension = 3;
const TinyVector<Dimension, uint64_t> node_size = [&] {
@@ -586,6 +591,10 @@ 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];
};
+ std::cout << "LogicalConnectivityBuilder:\n - prepare t = " << t << '\n';
+ t.stop();
+ t.start();
+
descriptor.cell_to_node_vector.resize(number_of_cells);
constexpr size_t nb_node_per_cell = 1 << Dimension;
for (size_t j = 0; j < number_of_cells; ++j) {
@@ -603,14 +612,36 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
descriptor.cell_to_node_vector[j][7] = node_number(cell_index + TinyVector<Dimension, uint64_t>{0, 1, 1});
}
}
+ std::cout << " - descriptor.cell_to_node_vector t = " << t << '\n';
+ t.stop();
+ t.start();
ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+ std::cout << " - cell->face, face->node connectivities t = " << t << '\n';
+ t.stop();
+
+ t.start();
ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
+ std::cout << " - face->edge, edge->node, cell->edge connectivities t = " << t << '\n';
+ t.stop();
+
+ t.start();
this->_buildBoundaryNodeList(cell_size, descriptor);
+ std::cout << " - boundaryNodeList t = " << t << '\n';
+ t.stop();
+
+ t.start();
this->_buildBoundaryEdgeList(cell_size, descriptor);
+ std::cout << " - boundaryEdgeList t = " << t << '\n';
+ t.stop();
+
+ t.start();
this->_buildBoundaryFaceList(cell_size, descriptor);
+ std::cout << " - boundaryFaceList t = " << t << '\n';
+ t.stop();
+ t.start();
descriptor.cell_owner_vector.resize(number_of_cells);
std::fill(descriptor.cell_owner_vector.begin(), descriptor.cell_owner_vector.end(), parallel::rank());
@@ -622,8 +653,12 @@ LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& ce
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ std::cout << " - final infos t = " << t << '\n';
+ t.stop();
+ t.start();
m_connectivity = Connectivity<Dimension>::build(descriptor);
+ std::cout << " - Connectivity<Dimension>::build(descriptor) t = " << t << '\n';
}
template <size_t Dimension>
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());