From c28b794bd647022b6be80992d4a0296fd9fd9bdc Mon Sep 17 00:00:00 2001
From: Stephane Del Pino <stephane.delpino44@gmail.com>
Date: Wed, 1 Jul 2020 19:20:23 +0200
Subject: [PATCH] Separate connectivity and mesh builders
Actually, classes are created but files remain to be split
---
src/language/modules/MeshModule.cpp | 2 +-
src/mesh/CMakeLists.txt | 2 +-
src/mesh/CartesianMeshBuilder.cpp | 290 ++--
src/mesh/CartesianMeshBuilder.hpp | 32 +-
...cpp => DiamondDualConnectivityBuilder.cpp} | 141 +-
src/mesh/DiamondDualConnectivityBuilder.hpp | 43 +
src/mesh/DiamondDualMeshBuilder.hpp | 28 -
src/mesh/GmshReader.cpp | 1469 +++++++++--------
src/mesh/GmshReader.hpp | 138 +-
src/mesh/MeshBuilderBase.cpp | 20 +-
src/mesh/MeshBuilderBase.hpp | 43 +-
11 files changed, 1286 insertions(+), 922 deletions(-)
rename src/mesh/{DiamondDualMeshBuilder.cpp => DiamondDualConnectivityBuilder.cpp} (83%)
create mode 100644 src/mesh/DiamondDualConnectivityBuilder.hpp
delete mode 100644 src/mesh/DiamondDualMeshBuilder.hpp
diff --git a/src/language/modules/MeshModule.cpp b/src/language/modules/MeshModule.cpp
index b9a657296..0181e5b2c 100644
--- a/src/language/modules/MeshModule.cpp
+++ b/src/language/modules/MeshModule.cpp
@@ -9,7 +9,7 @@
#include <language/utils/TypeDescriptor.hpp>
#include <mesh/CartesianMeshBuilder.hpp>
#include <mesh/Connectivity.hpp>
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
#include <mesh/GmshReader.hpp>
#include <mesh/Mesh.hpp>
#include <utils/Exceptions.hpp>
diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index 46a16a505..97b166608 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -3,7 +3,7 @@
add_library(
PugsMesh
CartesianMeshBuilder.cpp
- DiamondDualMeshBuilder.cpp
+ DiamondDualConnectivityBuilder.cpp
Connectivity.cpp
ConnectivityComputer.cpp
ConnectivityDispatcher.cpp
diff --git a/src/mesh/CartesianMeshBuilder.cpp b/src/mesh/CartesianMeshBuilder.cpp
index a272dde4d..caff700fa 100644
--- a/src/mesh/CartesianMeshBuilder.cpp
+++ b/src/mesh/CartesianMeshBuilder.cpp
@@ -7,16 +7,171 @@
#include <utils/Array.hpp>
#include <utils/Messenger.hpp>
+template <size_t Dimension>
+NodeValue<TinyVector<Dimension>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<Dimension>&,
+ const TinyVector<Dimension>&,
+ const TinyVector<Dimension, uint64_t>&,
+ const IConnectivity&) const
+{
+ static_assert(Dimension <= 3, "invalid dimension");
+}
+
+template <>
+NodeValue<TinyVector<1>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<1>& a,
+ const TinyVector<1>& b,
+ const TinyVector<1, uint64_t>& cell_size,
+ const IConnectivity& connectivity) const
+{
+ const double h = (b[0] - a[0]) / cell_size[0];
+ NodeValue<TinyVector<1>> xr(connectivity);
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) { xr[r] = a + r * h; });
+
+ return xr;
+}
+
+template <>
+NodeValue<TinyVector<2>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<2>& a,
+ const TinyVector<2>& b,
+ const TinyVector<2, uint64_t>& cell_size,
+ const IConnectivity& connectivity) const
+{
+ const TinyVector<2> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1]};
+
+ const TinyVector<2, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1};
+
+ const auto node_logic_id = [&](size_t r) {
+ const uint64_t r0 = r / node_size[1];
+ const uint64_t r1 = r % node_size[1];
+ return TinyVector<2, uint64_t>{r0, r1};
+ };
+
+ NodeValue<TinyVector<2>> xr(connectivity);
+
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const TinyVector<2, uint64_t> node_index = node_logic_id(r);
+ for (size_t i = 0; i < 2; ++i) {
+ xr[r][i] = a[i] + node_index[i] * h[i];
+ }
+ });
+
+ return xr;
+}
+
+template <>
+NodeValue<TinyVector<3>>
+CartesianMeshBuilder::_getNodeCoordinates(const TinyVector<3>& a,
+ const TinyVector<3>& b,
+ const TinyVector<3, uint64_t>& cell_size,
+ const IConnectivity& connectivity) const
+{
+ const TinyVector<3, uint64_t> node_size = [&] {
+ TinyVector node_size{cell_size};
+ for (size_t i = 0; i < 3; ++i) {
+ node_size[i] += 1;
+ }
+ return node_size;
+ }();
+
+ const auto node_logic_id = [&](size_t r) {
+ const size_t slice1 = node_size[1] * node_size[2];
+ const size_t& slice2 = node_size[2];
+ const uint64_t r0 = r / slice1;
+ const uint64_t r1 = (r - r0 * slice1) / slice2;
+ const uint64_t r2 = r - (r0 * slice1 + r1 * slice2);
+ return TinyVector<3, uint64_t>{r0, r1, r2};
+ };
+
+ const TinyVector<3> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1], (b[2] - a[2]) / cell_size[2]};
+
+ NodeValue<TinyVector<3>> xr(connectivity);
+ parallel_for(
+ connectivity.numberOfNodes(), PUGS_LAMBDA(NodeId r) {
+ const TinyVector<3, uint64_t> node_index = node_logic_id(r);
+ for (size_t i = 0; i < 3; ++i) {
+ xr[r][i] = a[i] + node_index[i] * h[i];
+ }
+ });
+
+ return xr;
+}
+
+template <size_t Dimension>
+void
+CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& cell_size)
+{
+ static_assert(Dimension <= 3, "unexpected dimension");
+
+ LogicalConnectivityBuilder logical_connectivity_builder{cell_size};
+
+ using ConnectivityType = Connectivity<Dimension>;
+
+ std::shared_ptr<const ConnectivityType> p_connectivity =
+ std::dynamic_pointer_cast<const ConnectivityType>(logical_connectivity_builder.connectivity());
+ const ConnectivityType& connectivity = *p_connectivity;
+
+ NodeValue<TinyVector<Dimension>> xr = _getNodeCoordinates(a, b, cell_size, connectivity);
+
+ m_mesh = std::make_shared<Mesh<ConnectivityType>>(p_connectivity, xr);
+}
+
+template <size_t Dimension>
+CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size)
+{
+ if (parallel::rank() == 0) {
+ TinyVector lenght = b - a;
+ for (size_t i = 0; i < Dimension; ++i) {
+ if (lenght[i] == 0) {
+ throw NormalError("invalid box definition corners share a component");
+ }
+ }
+
+ TinyVector<Dimension> corner0 = a;
+ TinyVector<Dimension> corner1 = b;
+
+ for (size_t i = 0; i < Dimension; ++i) {
+ if (corner0[i] > corner1[i]) {
+ std::swap(corner0[i], corner1[i]);
+ }
+ }
+
+ this->_buildCartesianMesh(corner0, corner1, size);
+ }
+ this->_dispatch<Dimension>();
+}
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
+ const TinyVector<1>&,
+ const TinyVector<1, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
+ const TinyVector<2>&,
+ const TinyVector<2, uint64_t>&);
+
+template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
+ const TinyVector<3>&,
+ const TinyVector<3, uint64_t>&);
+
+/// =====================
+
template <size_t Dimension>
inline void
-CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
{
static_assert(Dimension <= 3, "unexpected dimension");
}
template <>
inline void
-CartesianMeshBuilder::_buildBoundaryNodeList(
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
const TinyVector<1, uint64_t>& cell_size, // clazy:exclude=function-args-by-value
ConnectivityDescriptor& descriptor)
{
@@ -35,7 +190,7 @@ CartesianMeshBuilder::_buildBoundaryNodeList(
template <>
void
-CartesianMeshBuilder::_buildBoundaryNodeList(
+LogicalConnectivityBuilder::_buildBoundaryNodeList(
const TinyVector<2, uint64_t>& cell_size, // clazy:exclude=function-args-by-value
ConnectivityDescriptor& descriptor)
{
@@ -72,8 +227,8 @@ CartesianMeshBuilder::_buildBoundaryNodeList(
template <>
void
-CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell_size,
- ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
{
const TinyVector<3, uint64_t> node_size{cell_size[0] + 1, cell_size[1] + 1, cell_size[2] + 1};
@@ -132,15 +287,15 @@ CartesianMeshBuilder::_buildBoundaryNodeList(const TinyVector<3, uint64_t>& cell
template <size_t Dimension>
void
-CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
{
static_assert(Dimension == 3, "unexpected dimension");
}
template <>
void
-CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
- ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
{
using Edge = ConnectivityFace<2>;
const auto& node_number_vector = descriptor.node_number_vector;
@@ -242,14 +397,14 @@ CartesianMeshBuilder::_buildBoundaryEdgeList(const TinyVector<3, uint64_t>& cell
template <size_t Dimension>
void
-CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>&, ConnectivityDescriptor&)
{
static_assert(Dimension >= 2 and Dimension <= 3, "unexpected dimension");
}
template <>
void
-CartesianMeshBuilder::_buildBoundaryFaceList(
+LogicalConnectivityBuilder::_buildBoundaryFaceList(
const TinyVector<2, uint64_t>& cell_size, // clazy:exclude=function-args-by-value
ConnectivityDescriptor& descriptor)
{
@@ -316,8 +471,8 @@ CartesianMeshBuilder::_buildBoundaryFaceList(
template <>
void
-CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
- ConnectivityDescriptor& descriptor)
+LogicalConnectivityBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell_size,
+ ConnectivityDescriptor& descriptor)
{
using Face = ConnectivityFace<3>;
@@ -425,9 +580,7 @@ CartesianMeshBuilder::_buildBoundaryFaceList(const TinyVector<3, uint64_t>& cell
template <>
void
-CartesianMeshBuilder::_buildCartesianMesh(
- const TinyVector<1>& a, // clazy:exclude=function-args-by-value
- const TinyVector<1>& b, // clazy:exclude=function-args-by-value
+LogicalConnectivityBuilder::_buildConnectivity(
const TinyVector<1, uint64_t>& cell_size) // clazy:exclude=function-args-by-value
{
const size_t number_of_cells = cell_size[0];
@@ -465,24 +618,12 @@ CartesianMeshBuilder::_buildCartesianMesh(
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
- std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
- Connectivity1D& connectivity = *p_connectivity;
-
- const double h = (b[0] - a[0]) / cell_size[0];
-
- NodeValue<TinyVector<1>> xr(connectivity);
- for (NodeId r = 0; r < number_of_nodes; ++r) {
- xr[r] = a + r * h;
- }
-
- m_mesh = std::make_shared<Mesh<Connectivity1D>>(p_connectivity, xr);
+ m_connectivity = Connectivity1D::build(descriptor);
}
template <>
void
-CartesianMeshBuilder::_buildCartesianMesh(
- const TinyVector<2>& a, // clazy:exclude=function-args-by-value
- const TinyVector<2>& b, // clazy:exclude=function-args-by-value
+LogicalConnectivityBuilder::_buildConnectivity(
const TinyVector<2, uint64_t>& cell_size) // clazy:exclude=function-args-by-value
{
constexpr size_t Dimension = 2;
@@ -506,12 +647,6 @@ CartesianMeshBuilder::_buildCartesianMesh(
descriptor.cell_type_vector.resize(number_of_cells);
std::fill(descriptor.cell_type_vector.begin(), descriptor.cell_type_vector.end(), CellType::Quadrangle);
- const auto node_logic_id = [&](size_t r) {
- const uint64_t r0 = r / node_size[1];
- const uint64_t r1 = r % node_size[1];
- return TinyVector<Dimension, uint64_t>{r0, r1};
- };
-
const auto node_number = [&](const TinyVector<Dimension, uint64_t> node_logic_id) {
return node_logic_id[0] * node_size[1] + node_logic_id[1];
};
@@ -535,7 +670,7 @@ CartesianMeshBuilder::_buildCartesianMesh(
}
}
- MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
this->_buildBoundaryNodeList(cell_size, descriptor);
this->_buildBoundaryFaceList(cell_size, descriptor);
@@ -549,27 +684,12 @@ CartesianMeshBuilder::_buildCartesianMesh(
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
- std::shared_ptr p_connectivity = Connectivity<Dimension>::build(descriptor);
- Connectivity<Dimension>& connectivity = *p_connectivity;
-
- const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1]};
-
- NodeValue<TinyVector<Dimension>> xr(connectivity);
- for (NodeId r = 0; r < number_of_nodes; ++r) {
- const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
- for (size_t i = 0; i < Dimension; ++i) {
- xr[r][i] = a[i] + node_index[i] * h[i];
- }
- }
-
- m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+ m_connectivity = Connectivity<Dimension>::build(descriptor);
}
template <>
void
-CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
- const TinyVector<3>& b,
- const TinyVector<3, uint64_t>& cell_size)
+LogicalConnectivityBuilder::_buildConnectivity(const TinyVector<3, uint64_t>& cell_size)
{
constexpr size_t Dimension = 3;
@@ -618,15 +738,6 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
return TinyVector<Dimension, uint64_t>{j0, j1, j2};
};
- const auto node_logic_id = [&](size_t r) {
- const size_t slice1 = node_size[1] * node_size[2];
- const size_t& slice2 = node_size[2];
- const uint64_t r0 = r / slice1;
- const uint64_t r1 = (r - r0 * slice1) / slice2;
- const uint64_t r2 = r - (r0 * slice1 + r1 * slice2);
- return TinyVector<Dimension, uint64_t>{r0, r1, r2};
- };
-
const auto node_number = [&](const TinyVector<Dimension, uint64_t>& node_logic_id) {
return (node_logic_id[0] * node_size[1] + node_logic_id[1]) * node_size[2] + node_logic_id[2];
};
@@ -649,8 +760,8 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
}
}
- MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
- MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
+ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(descriptor);
+ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(descriptor);
this->_buildBoundaryNodeList(cell_size, descriptor);
this->_buildBoundaryEdgeList(cell_size, descriptor);
@@ -668,58 +779,19 @@ CartesianMeshBuilder::_buildCartesianMesh(const TinyVector<3>& a,
descriptor.node_owner_vector.resize(descriptor.node_number_vector.size());
std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
- std::shared_ptr p_connectivity = Connectivity<Dimension>::build(descriptor);
- Connectivity<Dimension>& connectivity = *p_connectivity;
-
- const TinyVector<Dimension> h{(b[0] - a[0]) / cell_size[0], (b[1] - a[1]) / cell_size[1],
- (b[2] - a[2]) / cell_size[2]};
-
- NodeValue<TinyVector<Dimension>> xr(connectivity);
- for (NodeId r = 0; r < number_of_nodes; ++r) {
- const TinyVector<Dimension, uint64_t> node_index = node_logic_id(r);
- for (size_t i = 0; i < Dimension; ++i) {
- xr[r][i] = a[i] + node_index[i] * h[i];
- }
- }
-
- m_mesh = std::make_shared<Mesh<Connectivity<Dimension>>>(p_connectivity, xr);
+ m_connectivity = Connectivity<Dimension>::build(descriptor);
}
template <size_t Dimension>
-CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<Dimension>& a,
- const TinyVector<Dimension>& b,
- const TinyVector<Dimension, uint64_t>& size)
+LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size)
{
if (parallel::rank() == 0) {
- TinyVector lenght = b - a;
- for (size_t i = 0; i < Dimension; ++i) {
- if (lenght[i] == 0) {
- throw NormalError("invalid box definition corners share a component");
- }
- }
-
- TinyVector<Dimension> corner0 = a;
- TinyVector<Dimension> corner1 = b;
-
- for (size_t i = 0; i < Dimension; ++i) {
- if (corner0[i] > corner1[i]) {
- std::swap(corner0[i], corner1[i]);
- }
- }
-
- this->_buildCartesianMesh(corner0, corner1, size);
+ this->_buildConnectivity(size);
}
- this->_dispatch<Dimension>();
}
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<1>&,
- const TinyVector<1>&,
- const TinyVector<1, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<1, uint64_t>&);
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<2>&,
- const TinyVector<2>&,
- const TinyVector<2, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<2, uint64_t>&);
-template CartesianMeshBuilder::CartesianMeshBuilder(const TinyVector<3>&,
- const TinyVector<3>&,
- const TinyVector<3, uint64_t>&);
+template LogicalConnectivityBuilder::LogicalConnectivityBuilder(const TinyVector<3, uint64_t>&);
diff --git a/src/mesh/CartesianMeshBuilder.hpp b/src/mesh/CartesianMeshBuilder.hpp
index 712511583..9ed1ea81c 100644
--- a/src/mesh/CartesianMeshBuilder.hpp
+++ b/src/mesh/CartesianMeshBuilder.hpp
@@ -8,6 +8,28 @@
#include <memory>
class CartesianMeshBuilder : public MeshBuilderBase
+{
+ private:
+ template <size_t Dimension>
+ void _buildCartesianMesh(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size);
+
+ template <size_t Dimension>
+ NodeValue<TinyVector<Dimension>> _getNodeCoordinates(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size,
+ const IConnectivity& connectivity) const;
+
+ public:
+ template <size_t Dimension>
+ CartesianMeshBuilder(const TinyVector<Dimension>& a,
+ const TinyVector<Dimension>& b,
+ const TinyVector<Dimension, uint64_t>& size);
+ ~CartesianMeshBuilder() = default;
+};
+
+class LogicalConnectivityBuilder : public ConnectivityBuilderBase
{
private:
template <size_t Dimension>
@@ -20,16 +42,12 @@ class CartesianMeshBuilder : public MeshBuilderBase
void _buildBoundaryFaceList(const TinyVector<Dimension, uint64_t>& cell_size, ConnectivityDescriptor& descriptor);
template <size_t Dimension>
- void _buildCartesianMesh(const TinyVector<Dimension>& a,
- const TinyVector<Dimension>& b,
- const TinyVector<Dimension, uint64_t>& size);
+ void _buildConnectivity(const TinyVector<Dimension, uint64_t>& size);
public:
template <size_t Dimension>
- CartesianMeshBuilder(const TinyVector<Dimension>& a,
- const TinyVector<Dimension>& b,
- const TinyVector<Dimension, uint64_t>& size);
- ~CartesianMeshBuilder() = default;
+ LogicalConnectivityBuilder(const TinyVector<Dimension, uint64_t>& size);
+ ~LogicalConnectivityBuilder() = default;
};
#endif // CARTESIAN_MESH_BUILDER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.cpp b/src/mesh/DiamondDualConnectivityBuilder.cpp
similarity index 83%
rename from src/mesh/DiamondDualMeshBuilder.cpp
rename to src/mesh/DiamondDualConnectivityBuilder.cpp
index 4c8c06580..585982245 100644
--- a/src/mesh/DiamondDualMeshBuilder.cpp
+++ b/src/mesh/DiamondDualConnectivityBuilder.cpp
@@ -1,4 +1,4 @@
-#include <mesh/DiamondDualMeshBuilder.hpp>
+#include <mesh/DiamondDualConnectivityBuilder.hpp>
#include <mesh/Connectivity.hpp>
#include <mesh/ConnectivityDescriptor.hpp>
@@ -12,8 +12,8 @@
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
- ConnectivityDescriptor& diamond_descriptor)
+DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<Dimension>& primal_connectivity,
+ ConnectivityDescriptor& diamond_descriptor)
{
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
@@ -154,8 +154,8 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor(const Connectivity<D
template <>
void
-DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
- ConnectivityDescriptor& diamond_descriptor)
+DiamondDualConnectivityBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivity<1>& primal_connectivity,
+ ConnectivityDescriptor& diamond_descriptor)
{
const size_t primal_number_of_nodes = primal_connectivity.numberOfNodes();
const size_t primal_number_of_faces = primal_connectivity.numberOfFaces();
@@ -227,7 +227,7 @@ DiamondDualMeshBuilder::_buildDiamondConnectivityDescriptor<1>(const Connectivit
template <size_t Dimension>
void
-DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
+DiamondDualConnectivityBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>& p_mesh)
{
static_assert(Dimension <= 3, "invalid mesh dimension");
using ConnectivityType = Connectivity<Dimension>;
@@ -242,9 +242,9 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
this->_buildDiamondConnectivityDescriptor(primal_connectivity, diamond_descriptor);
if constexpr (Dimension > 1) {
- MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
+ ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<Dimension>(diamond_descriptor);
if constexpr (Dimension > 2) {
- MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
+ ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<Dimension>(diamond_descriptor);
}
}
@@ -474,57 +474,128 @@ DiamondDualMeshBuilder::_buildDiamondMeshFrom(const std::shared_ptr<const IMesh>
}
}
- std::shared_ptr p_diamond_connectivity = ConnectivityType::build(diamond_descriptor);
- ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+ m_connectivity = ConnectivityType::build(diamond_descriptor);
+}
+
+DiamondDualConnectivityBuilder::DiamondDualConnectivityBuilder(const std::shared_ptr<const IMesh>& p_mesh)
+{
+ switch (p_mesh->dimension()) {
+ case 1: {
+ this->_buildDiamondMeshFrom<1>(p_mesh);
+ break;
+ }
+ case 2: {
+ this->_buildDiamondMeshFrom<2>(p_mesh);
+ break;
+ }
+ case 3: {
+ this->_buildDiamondMeshFrom<3>(p_mesh);
+ break;
+ }
+ default: {
+ throw UnexpectedError("invalid mesh dimension: " + std::to_string(p_mesh->dimension()));
+ }
+ }
+}
+
+template <size_t Dimension>
+void
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(
+ const Mesh<Connectivity<Dimension>>& primal_mesh,
+ const std::shared_ptr<const Connectivity<Dimension>>& p_diamond_connectivity)
+{
+ using ConnectivityType = Connectivity<Dimension>;
+ using MeshType = Mesh<Connectivity<Dimension>>;
+
+ const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
NodeValue<TinyVector<Dimension>> diamond_xr{diamond_connectivity};
- const auto primal_xr = primal_mesh.xr();
+ const NodeValue<const TinyVector<Dimension>> primal_xr = primal_mesh.xr();
MeshData<MeshType> primal_mesh_data{primal_mesh};
- const auto primal_xj = primal_mesh_data.xj();
+ const CellValue<const TinyVector<Dimension>> primal_xj = primal_mesh_data.xj();
- {
- if constexpr (Dimension == 1) {
- const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
- NodeId next_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) {
- diamond_xr[next_node_id++] = primal_xr[primal_node_id];
- }
- }
+ NodeId i_node = 0;
+ for (; i_node < primal_mesh.numberOfNodes(); ++i_node) {
+ diamond_xr[i_node] = primal_xr[i_node];
+ }
- for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
- diamond_xr[next_node_id++] = primal_xj[i_cell];
- }
+ for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
+ diamond_xr[i_node++] = primal_xj[i_cell];
+ }
- } else {
- NodeId i_node = 0;
- for (; i_node < primal_number_of_nodes; ++i_node) {
- diamond_xr[i_node] = primal_xr[i_node];
- }
+ m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
+}
- for (CellId i_cell = 0; i_cell < primal_number_of_cells; ++i_cell) {
- diamond_xr[i_node++] = primal_xj[i_cell];
- }
+template <>
+void
+DiamondDualMeshBuilder::_buildDualDiamondMeshFrom(const Mesh<Connectivity<1>>& primal_mesh,
+ const std::shared_ptr<const Connectivity<1>>& p_diamond_connectivity)
+{
+ using ConnectivityType = Connectivity<1>;
+ using MeshType = Mesh<Connectivity<1>>;
+
+ const ConnectivityType& primal_connectivity = primal_mesh.connectivity();
+ const auto& node_to_cell_matrix = primal_connectivity.nodeToCellMatrix();
+
+ const ConnectivityType& diamond_connectivity = *p_diamond_connectivity;
+
+ NodeValue<TinyVector<1>> diamond_xr{diamond_connectivity};
+
+ const NodeValue<const TinyVector<1>> primal_xr = primal_mesh.xr();
+ MeshData<MeshType> primal_mesh_data{primal_mesh};
+ const CellValue<const TinyVector<1>> primal_xj = primal_mesh_data.xj();
+
+ NodeId next_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) {
+ diamond_xr[next_node_id++] = primal_xr[primal_node_id];
}
}
+ for (CellId i_cell = 0; i_cell < primal_mesh.numberOfCells(); ++i_cell) {
+ diamond_xr[next_node_id++] = primal_xj[i_cell];
+ }
+
m_mesh = std::make_shared<MeshType>(p_diamond_connectivity, diamond_xr);
}
DiamondDualMeshBuilder::DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>& p_mesh)
{
+ DiamondDualConnectivityBuilder builder(p_mesh);
+
switch (p_mesh->dimension()) {
case 1: {
- this->_buildDiamondMeshFrom<1>(p_mesh);
+ using ConnectivityType = Connectivity<1>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+ this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
break;
}
case 2: {
- this->_buildDiamondMeshFrom<2>(p_mesh);
+ using ConnectivityType = Connectivity<2>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+ this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
break;
}
case 3: {
- this->_buildDiamondMeshFrom<3>(p_mesh);
+ using ConnectivityType = Connectivity<3>;
+ using MeshType = Mesh<ConnectivityType>;
+
+ std::shared_ptr p_connectivity = std::dynamic_pointer_cast<const ConnectivityType>(builder.connectivity());
+
+ const MeshType& mesh = dynamic_cast<const MeshType&>(*p_mesh);
+
+ this->_buildDualDiamondMeshFrom(mesh, p_connectivity);
break;
}
default: {
diff --git a/src/mesh/DiamondDualConnectivityBuilder.hpp b/src/mesh/DiamondDualConnectivityBuilder.hpp
new file mode 100644
index 000000000..a2fe946b9
--- /dev/null
+++ b/src/mesh/DiamondDualConnectivityBuilder.hpp
@@ -0,0 +1,43 @@
+#ifndef DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
+#define DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
+
+#include <mesh/IMesh.hpp>
+#include <mesh/MeshBuilderBase.hpp>
+
+#include <memory>
+
+template <size_t>
+class Connectivity;
+
+template <typename ConnectivityType>
+class Mesh;
+
+class ConnectivityDescriptor;
+
+class DiamondDualConnectivityBuilder : public ConnectivityBuilderBase
+{
+ private:
+ template <size_t Dimension>
+ void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
+
+ template <size_t Dimension>
+ void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
+
+ public:
+ DiamondDualConnectivityBuilder(const std::shared_ptr<const IMesh>&);
+ ~DiamondDualConnectivityBuilder() = default;
+};
+
+class DiamondDualMeshBuilder : public MeshBuilderBase
+{
+ private:
+ template <size_t Dimension>
+ void _buildDualDiamondMeshFrom(const Mesh<Connectivity<Dimension>>&,
+ const std::shared_ptr<const Connectivity<Dimension>>&);
+
+ public:
+ DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
+ ~DiamondDualMeshBuilder() = default;
+};
+
+#endif // DIAMOND_DUAL_CONNECTIVITY_BUILDER_HPP
diff --git a/src/mesh/DiamondDualMeshBuilder.hpp b/src/mesh/DiamondDualMeshBuilder.hpp
deleted file mode 100644
index a6af1b830..000000000
--- a/src/mesh/DiamondDualMeshBuilder.hpp
+++ /dev/null
@@ -1,28 +0,0 @@
-#ifndef DIAMOND_DUAL_MESH_BUILDER_HPP
-#define DIAMOND_DUAL_MESH_BUILDER_HPP
-
-#include <mesh/IMesh.hpp>
-#include <mesh/MeshBuilderBase.hpp>
-
-#include <memory>
-
-template <size_t>
-class Connectivity;
-
-class ConnectivityDescriptor;
-
-class DiamondDualMeshBuilder : public MeshBuilderBase
-{
- private:
- template <size_t Dimension>
- void _buildDiamondConnectivityDescriptor(const Connectivity<Dimension>&, ConnectivityDescriptor&);
-
- template <size_t Dimension>
- void _buildDiamondMeshFrom(const std::shared_ptr<const IMesh>&);
-
- public:
- DiamondDualMeshBuilder(const std::shared_ptr<const IMesh>&);
- ~DiamondDualMeshBuilder() = default;
-};
-
-#endif // DIAMOND_DUAL_MESH_BUILDER_HPP
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index e99c04722..6ca41953f 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -23,6 +23,450 @@
#include <sstream>
#include <unordered_map>
+template <size_t Dimension>
+class GmshConnectivityBuilder : public ConnectivityBuilderBase
+{
+ public:
+ GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells);
+};
+
+template <>
+GmshConnectivityBuilder<1>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+ 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);
+
+ 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];
+ }
+ descriptor.cell_type_vector[j] = CellType::Line;
+ descriptor.cell_number_vector[j] = gmsh_data.__edges_number[j];
+ }
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+ for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
+ const unsigned int point_number = gmsh_data.__points[r];
+ const unsigned int& ref = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+ m_connectivity = Connectivity1D::build(descriptor);
+}
+
+template <>
+GmshConnectivityBuilder<2>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+ 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);
+
+ 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];
+ }
+ descriptor.cell_type_vector[j] = CellType::Triangle;
+ descriptor.cell_number_vector[j] = gmsh_data.__triangles_number[j];
+ }
+
+ 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];
+ }
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+ for (unsigned int r = 0; r < gmsh_data.__triangles_ref.size(); ++r) {
+ const unsigned int elem_number = gmsh_data.__triangles_ref[r];
+ const unsigned int& ref = gmsh_data.__triangles_ref[r];
+ ref_cells_map[ref].push_back(elem_number);
+ }
+
+ for (unsigned int j = 0; j < gmsh_data.__quadrangles_ref.size(); ++j) {
+ const size_t elem_number = nb_triangles + j;
+ const unsigned int& ref = gmsh_data.__quadrangles_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_cell_list.first);
+ descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
+ }
+
+ 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 = [&] {
+ 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 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& 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 (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]);
+
+ 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[descriptor.face_number_vector[edge_id]] = edge_id;
+ face_number_id_map[descriptor.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;
+ }
+ }
+ }
+
+ 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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_data.m_physical_ref_map.at(ref_face_list.first);
+ descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
+ }
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
+ for (unsigned int r = 0; r < gmsh_data.__points.size(); ++r) {
+ const unsigned int point_number = gmsh_data.__points[r];
+ const unsigned int& ref = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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.node_owner_vector.resize(descriptor.node_number_vector.size());
+ std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+
+ m_connectivity = Connectivity2D::build(descriptor);
+}
+
+template <>
+GmshConnectivityBuilder<3>::GmshConnectivityBuilder(const GmshReader::GmshData& gmsh_data, const size_t nb_cells)
+{
+ 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);
+
+ 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];
+ }
+ 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_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];
+ }
+
+ std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
+ for (unsigned int r = 0; r < gmsh_data.__tetrahedra_ref.size(); ++r) {
+ const unsigned int elem_number = gmsh_data.__tetrahedra_ref[r];
+ const unsigned int& ref = gmsh_data.__tetrahedra_ref[r];
+ ref_cells_map[ref].push_back(elem_number);
+ }
+
+ for (unsigned int j = 0; j < gmsh_data.__hexahedra_ref.size(); ++j) {
+ const size_t elem_number = nb_tetrahedra + j;
+ const unsigned int& ref = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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 < 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& ref = gmsh_data.__triangles_ref[f];
+ ref_faces_map[ref].push_back(face_id);
+
+ if (descriptor.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]);
+
+ 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] = gmsh_data.__quadrangles_number[f];
+ face_number_id_map[descriptor.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& 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 (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]);
+
+ 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] = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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 < 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];
+ ref_edges_map[ref].push_back(edge_id);
+
+ if (descriptor.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]);
+
+ 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] = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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 < gmsh_data.__points.size(); ++r) {
+ const unsigned int point_number = gmsh_data.__points[r];
+ const unsigned int& ref = gmsh_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 GmshReader::PhysicalRefId& physical_ref_id = gmsh_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());
+
+ m_connectivity = Connectivity3D::build(descriptor);
+}
+
GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
{
if (parallel::rank() == 0) {
@@ -106,7 +550,6 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
throw NormalError("Cannot read Gmsh format '" + std::to_string(fileVersion) + "'");
}
int fileType = this->_getInteger();
- __binary = (fileType == 1);
if ((fileType < 0) or (fileType > 1)) {
throw NormalError("Cannot read Gmsh file type '" + std::to_string(fileType) + "'");
}
@@ -116,10 +559,6 @@ GmshReader::GmshReader(const std::string& filename) : m_filename(filename)
throw NormalError("Data size not supported '" + std::to_string(dataSize) + "'");
}
- if (__binary) {
- throw NotImplementedError("cannot read binary files");
- }
-
kw = this->__nextKeyword();
if (kw.second != ENDMESHFORMAT) {
throw NormalError("reading file '" + m_filename + "': expecting $EndMeshFormat, '" + kw.first + "' was found");
@@ -189,71 +628,18 @@ GmshReader::__readVertices()
throw NormalError("reading file '" + this->m_filename + "': number of vertices is negative");
}
- __verticesNumbers.resize(numberOfVerices);
- __vertices = Array<R3>(numberOfVerices);
+ m_mesh_data.__verticesNumbers.resize(numberOfVerices);
+ m_mesh_data.__vertices = Array<R3>(numberOfVerices);
- if (not __binary) {
- for (int i = 0; i < numberOfVerices; ++i) {
- __verticesNumbers[i] = this->_getInteger();
- const double x = this->_getReal();
- const double y = this->_getReal();
- const double z = this->_getReal();
- __vertices[i] = TinyVector<3, double>(x, y, z);
- }
- } else {
- throw NotImplementedError("cannot read binary format");
+ for (int i = 0; i < numberOfVerices; ++i) {
+ m_mesh_data.__verticesNumbers[i] = this->_getInteger();
+ const double x = this->_getReal();
+ const double y = this->_getReal();
+ const double z = this->_getReal();
+ m_mesh_data.__vertices[i] = TinyVector<3, double>(x, y, z);
}
}
-// void
-// GmshReader::__readElements1()
-// {
-// const int numberOfElements = this->_getInteger();
-// std::cout << "- Number Of Elements: " << numberOfElements << '\n';
-// if (numberOfElements<0) {
-// throw ErrorHandler(__FILE__,__LINE__,
-// "reading file '"+m_filename
-// +"': number of elements is negative",
-// ErrorHandler::normal);
-// }
-
-// __elementType.resize(numberOfElements);
-// __references.resize(numberOfElements);
-// __elementVertices.resize(numberOfElements);
-
-// for (int i=0; i<numberOfElements; ++i) {
-// this->_getInteger(); // drops element number
-// const int elementType = this->_getInteger()-1;
-
-// if ((elementType < 0) or (elementType > 14)) {
-// throw ErrorHandler(__FILE__,__LINE__,
-// "reading file '"+m_filename
-// +"': unknown element type
-// '"+std::to_string(elementType)+"'", ErrorHandler::normal);
-// }
-// __elementType[i] = elementType;
-// __references[i] = this->_getInteger(); // physical reference
-// this->_getInteger(); // drops entity reference
-
-// const int numberOfVertices = this->_getInteger();
-// if (numberOfVertices != __numberOfPrimitiveNodes[elementType]) {
-// std::stringstream errorMsg;
-// errorMsg << "reading file '" <<m_filename
-// << "':number of vertices (" << numberOfVertices
-// << ") does not match expected ("
-// << __numberOfPrimitiveNodes[elementType] << ")" << std::ends;
-
-// throw ErrorHandler(__FILE__,__LINE__,
-// errorMsg.str(),
-// ErrorHandler::normal);
-// }
-// __elementVertices[i].resize(numberOfVertices);
-// for (int j=0; j<numberOfVertices; ++j) {
-// __elementVertices[i][j] = this->_getInteger();
-// }
-// }
-// }
-
void
GmshReader::__readElements2_2()
{
@@ -263,96 +649,31 @@ GmshReader::__readElements2_2()
throw NormalError("reading file '" + m_filename + "': number of elements is negative");
}
- __elementNumber.resize(numberOfElements);
- __elementType.resize(numberOfElements);
- __references.resize(numberOfElements);
- __elementVertices.resize(numberOfElements);
-
- if (not __binary) {
- for (int i = 0; i < numberOfElements; ++i) {
- __elementNumber[i] = this->_getInteger();
- const int elementType = this->_getInteger() - 1;
+ m_mesh_data.__elementNumber.resize(numberOfElements);
+ m_mesh_data.__elementType.resize(numberOfElements);
+ m_mesh_data.__references.resize(numberOfElements);
+ m_mesh_data.__elementVertices.resize(numberOfElements);
- if ((elementType < 0) or (elementType > 14)) {
- throw NormalError("reading file '" + m_filename + "': unknown element type '" + std::to_string(elementType) +
- "'");
- }
- __elementType[i] = elementType;
- const int numberOfTags = this->_getInteger();
- __references[i] = this->_getInteger(); // physical reference
- for (int tag = 1; tag < numberOfTags; ++tag) {
- this->_getInteger(); // drops remaining tags
- }
+ for (int i = 0; i < numberOfElements; ++i) {
+ m_mesh_data.__elementNumber[i] = this->_getInteger();
+ const int elementType = this->_getInteger() - 1;
- const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
- __elementVertices[i].resize(numberOfVertices);
- for (int j = 0; j < numberOfVertices; ++j) {
- __elementVertices[i][j] = this->_getInteger();
- }
+ if ((elementType < 0) or (elementType > 14)) {
+ throw NormalError("reading file '" + m_filename + "': unknown element type '" + std::to_string(elementType) +
+ "'");
+ }
+ m_mesh_data.__elementType[i] = elementType;
+ const int numberOfTags = this->_getInteger();
+ m_mesh_data.__references[i] = this->_getInteger(); // physical reference
+ for (int tag = 1; tag < numberOfTags; ++tag) {
+ this->_getInteger(); // drops remaining tags
}
- } else {
- // fseek(__ifh,1L,SEEK_CUR);
- // int i=0;
- // for (;i<numberOfElements;) {
- // int elementType = 0;
- // fread(reinterpret_cast<char*>(&elementType),sizeof(int),1,__ifh);
- // if (__convertEndian) {
- // invertEndianess(elementType);
- // }
- // --elementType;
-
- // if ((elementType < 0) or (elementType > 14)) {
- // throw ErrorHandler(__FILE__,__LINE__,
- // "reading file '"+m_filename
- // +"': unknown element type
- // '"+std::to_string(elementType)+"'",
- // ErrorHandler::normal);
- // }
-
- // int elementTypeNumber = 0;
- // fread(reinterpret_cast<char*>(&elementTypeNumber),sizeof(int),1,__ifh);
- // if (__convertEndian) {
- // invertEndianess(elementTypeNumber);
- // }
-
- // int numberOfTags = 0;
- // fread(reinterpret_cast<char*>(&numberOfTags),sizeof(int),1,__ifh);
- // if (__convertEndian) {
- // invertEndianess(numberOfTags);
- // }
-
- // for (int k=0; k<elementTypeNumber; ++k) {
- // __elementType[i] = elementType;
-
- // int elementNumber = 0;
- // fread(reinterpret_cast<char*>(&elementNumber),sizeof(int),1,__ifh);
-
- // int reference = 0;
- // fread(reinterpret_cast<char*>(&reference),sizeof(int),1,__ifh);
-
- // __references[i] = reference; // physical reference
- // for (int tag=1; tag<numberOfTags; ++tag) {
- // int t;
- // fread(reinterpret_cast<char*>(&t),sizeof(int),1,__ifh);
- // }
-
- // const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
- // __elementVertices[i].resize(numberOfVertices);
- // fread(reinterpret_cast<char*>(&(__elementVertices[i][0])),sizeof(int),numberOfVertices,__ifh);
- // ++i;
- // }
- // }
- // if (__convertEndian) {
- // for (size_t i=0; i<__references.size(); ++i) {
- // invertEndianess(__references[i]);
- // }
- // for (size_t i=0; i<__elementVertices.size(); ++i) {
- // for (size_t j=0; j<__elementVertices[i].size(); ++i) {
- // invertEndianess(__elementVertices[i][j]);
- // }
- // }
- // }
+ const int numberOfVertices = __numberOfPrimitiveNodes[elementType];
+ m_mesh_data.__elementVertices[i].resize(numberOfVertices);
+ for (int j = 0; j < numberOfVertices; ++j) {
+ m_mesh_data.__elementVertices[i][j] = this->_getInteger();
+ }
}
}
@@ -369,27 +690,277 @@ GmshReader::__readPhysicalNames2_2()
PhysicalRefId physical_ref_id(physical_dimension, RefId(physical_number, physical_name));
- if (auto i_searched_physical_ref_id = m_physical_ref_map.find(physical_number);
- i_searched_physical_ref_id != m_physical_ref_map.end()) {
+ if (auto i_searched_physical_ref_id = m_mesh_data.m_physical_ref_map.find(physical_number);
+ i_searched_physical_ref_id != m_mesh_data.m_physical_ref_map.end()) {
std::stringstream os;
os << "Physical reference id '" << physical_ref_id << "' already defined as '"
<< i_searched_physical_ref_id->second << "'!";
throw NormalError(os.str());
}
- m_physical_ref_map[physical_number] = physical_ref_id;
+ m_mesh_data.m_physical_ref_map[physical_number] = physical_ref_id;
}
}
+// 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()
{
TinyVector<15, int> elementNumber = zero;
std::vector<std::set<size_t>> elementReferences(15);
- for (size_t i = 0; i < __elementType.size(); ++i) {
- const int elementType = __elementType[i];
+ for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
+ const int elementType = m_mesh_data.__elementType[i];
elementNumber[elementType]++;
- elementReferences[elementType].insert(__references[i]);
+ elementReferences[elementType].insert(m_mesh_data.__references[i]);
}
for (size_t i = 0; i < elementNumber.dimension(); ++i) {
@@ -409,39 +980,39 @@ GmshReader::__proceedData()
switch (i) {
// Supported entities
case 0: { // edges
- __edges = Array<Edge>(elementNumber[i]);
- __edges_ref.resize(elementNumber[i]);
- __edges_number.resize(elementNumber[i]);
+ m_mesh_data.__edges = Array<GmshData::Edge>(elementNumber[i]);
+ m_mesh_data.__edges_ref.resize(elementNumber[i]);
+ m_mesh_data.__edges_number.resize(elementNumber[i]);
break;
}
case 1: { // triangles
- __triangles = Array<Triangle>(elementNumber[i]);
- __triangles_ref.resize(elementNumber[i]);
- __triangles_number.resize(elementNumber[i]);
+ m_mesh_data.__triangles = Array<GmshData::Triangle>(elementNumber[i]);
+ m_mesh_data.__triangles_ref.resize(elementNumber[i]);
+ m_mesh_data.__triangles_number.resize(elementNumber[i]);
break;
}
case 2: { // quadrangles
- __quadrangles = Array<Quadrangle>(elementNumber[i]);
- __quadrangles_ref.resize(elementNumber[i]);
- __quadrangles_number.resize(elementNumber[i]);
+ m_mesh_data.__quadrangles = Array<GmshData::Quadrangle>(elementNumber[i]);
+ m_mesh_data.__quadrangles_ref.resize(elementNumber[i]);
+ m_mesh_data.__quadrangles_number.resize(elementNumber[i]);
break;
}
case 3: { // tetrahedra
- __tetrahedra = Array<Tetrahedron>(elementNumber[i]);
- __tetrahedra_ref.resize(elementNumber[i]);
- __tetrahedra_number.resize(elementNumber[i]);
+ m_mesh_data.__tetrahedra = Array<GmshData::Tetrahedron>(elementNumber[i]);
+ m_mesh_data.__tetrahedra_ref.resize(elementNumber[i]);
+ m_mesh_data.__tetrahedra_number.resize(elementNumber[i]);
break;
}
case 4: { // hexahedra
- __hexahedra = Array<Hexahedron>(elementNumber[i]);
- __hexahedra_ref.resize(elementNumber[i]);
- __hexahedra_number.resize(elementNumber[i]);
+ m_mesh_data.__hexahedra = Array<GmshData::Hexahedron>(elementNumber[i]);
+ m_mesh_data.__hexahedra_ref.resize(elementNumber[i]);
+ m_mesh_data.__hexahedra_number.resize(elementNumber[i]);
break;
}
case 14: { // point
- __points = Array<Point>(elementNumber[i]);
- __points_ref.resize(elementNumber[i]);
- __points_number.resize(elementNumber[i]);
+ m_mesh_data.__points = Array<GmshData::Point>(elementNumber[i]);
+ m_mesh_data.__points_ref.resize(elementNumber[i]);
+ m_mesh_data.__points_number.resize(elementNumber[i]);
break;
}
// Unsupported entities
@@ -466,8 +1037,8 @@ GmshReader::__proceedData()
std::cout << "- Building correspondance table\n";
int minNumber = std::numeric_limits<int>::max();
int maxNumber = std::numeric_limits<int>::min();
- for (size_t i = 0; i < __verticesNumbers.size(); ++i) {
- const int& vertexNumber = __verticesNumbers[i];
+ for (size_t i = 0; i < m_mesh_data.__verticesNumbers.size(); ++i) {
+ const int& vertexNumber = m_mesh_data.__verticesNumbers[i];
minNumber = std::min(minNumber, vertexNumber);
maxNumber = std::max(maxNumber, vertexNumber);
}
@@ -477,112 +1048,112 @@ GmshReader::__proceedData()
}
// A value of -1 means that the vertex is unknown
- __verticesCorrepondance.resize(maxNumber + 1, -1);
+ m_mesh_data.__verticesCorrepondance.resize(maxNumber + 1, -1);
- for (size_t i = 0; i < __verticesNumbers.size(); ++i) {
- __verticesCorrepondance[__verticesNumbers[i]] = i;
+ for (size_t i = 0; i < m_mesh_data.__verticesNumbers.size(); ++i) {
+ m_mesh_data.__verticesCorrepondance[m_mesh_data.__verticesNumbers[i]] = i;
}
// reset element number to count them while filling structures
elementNumber = zero;
// Element structures filling
- for (size_t i = 0; i < __elementType.size(); ++i) {
- std::vector<int>& elementVertices = __elementVertices[i];
- switch (__elementType[i]) {
+ for (size_t i = 0; i < m_mesh_data.__elementType.size(); ++i) {
+ std::vector<int>& elementVertices = m_mesh_data.__elementVertices[i];
+ switch (m_mesh_data.__elementType[i]) {
// Supported entities
case 0: { // edge
int& edgeNumber = elementNumber[0];
- const int a = __verticesCorrepondance[elementVertices[0]];
- const int b = __verticesCorrepondance[elementVertices[1]];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
if ((a < 0) or (b < 0)) {
throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
" [bad vertices definition]");
}
- __edges[edgeNumber] = Edge(a, b);
- __edges_ref[edgeNumber] = __references[i];
- __edges_number[edgeNumber] = __elementNumber[i];
+ m_mesh_data.__edges[edgeNumber] = GmshData::Edge(a, b);
+ m_mesh_data.__edges_ref[edgeNumber] = m_mesh_data.__references[i];
+ m_mesh_data.__edges_number[edgeNumber] = m_mesh_data.__elementNumber[i];
edgeNumber++; // one more edge
break;
}
case 1: { // triangles
int& triangleNumber = elementNumber[1];
- const int a = __verticesCorrepondance[elementVertices[0]];
- const int b = __verticesCorrepondance[elementVertices[1]];
- const int c = __verticesCorrepondance[elementVertices[2]];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+ const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
if ((a < 0) or (b < 0) or (c < 0)) {
throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
" [bad vertices definition]");
}
- __triangles[triangleNumber] = Triangle(a, b, c);
- __triangles_ref[triangleNumber] = __references[i];
- __triangles_number[triangleNumber] = __elementNumber[i];
+ m_mesh_data.__triangles[triangleNumber] = GmshData::Triangle(a, b, c);
+ m_mesh_data.__triangles_ref[triangleNumber] = m_mesh_data.__references[i];
+ m_mesh_data.__triangles_number[triangleNumber] = m_mesh_data.__elementNumber[i];
triangleNumber++; // one more triangle
break;
}
case 2: { // quadrangle
int& quadrilateralNumber = elementNumber[2];
- const int a = __verticesCorrepondance[elementVertices[0]];
- const int b = __verticesCorrepondance[elementVertices[1]];
- const int c = __verticesCorrepondance[elementVertices[2]];
- const int d = __verticesCorrepondance[elementVertices[3]];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+ const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+ const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
" [bad vertices definition]");
}
- __quadrangles[quadrilateralNumber] = Quadrangle(a, b, c, d);
- __quadrangles_ref[quadrilateralNumber] = __references[i];
- __quadrangles_number[quadrilateralNumber] = __elementNumber[i];
+ m_mesh_data.__quadrangles[quadrilateralNumber] = GmshData::Quadrangle(a, b, c, d);
+ m_mesh_data.__quadrangles_ref[quadrilateralNumber] = m_mesh_data.__references[i];
+ m_mesh_data.__quadrangles_number[quadrilateralNumber] = m_mesh_data.__elementNumber[i];
quadrilateralNumber++; // one more quadrangle
break;
}
case 3: { // tetrahedra
int& tetrahedronNumber = elementNumber[3];
- const int a = __verticesCorrepondance[elementVertices[0]];
- const int b = __verticesCorrepondance[elementVertices[1]];
- const int c = __verticesCorrepondance[elementVertices[2]];
- const int d = __verticesCorrepondance[elementVertices[3]];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+ const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+ const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
if ((a < 0) or (b < 0) or (c < 0) or (d < 0)) {
throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
" [bad vertices definition]");
}
- __tetrahedra[tetrahedronNumber] = Tetrahedron(a, b, c, d);
- __tetrahedra_ref[tetrahedronNumber] = __references[i];
- __tetrahedra_number[tetrahedronNumber] = __elementNumber[i];
+ m_mesh_data.__tetrahedra[tetrahedronNumber] = GmshData::Tetrahedron(a, b, c, d);
+ m_mesh_data.__tetrahedra_ref[tetrahedronNumber] = m_mesh_data.__references[i];
+ m_mesh_data.__tetrahedra_number[tetrahedronNumber] = m_mesh_data.__elementNumber[i];
tetrahedronNumber++; // one more tetrahedron
break;
}
case 4: { // hexaredron
int& hexahedronNumber = elementNumber[4];
- const int a = __verticesCorrepondance[elementVertices[0]];
- const int b = __verticesCorrepondance[elementVertices[1]];
- const int c = __verticesCorrepondance[elementVertices[2]];
- const int d = __verticesCorrepondance[elementVertices[3]];
- const int e = __verticesCorrepondance[elementVertices[4]];
- const int f = __verticesCorrepondance[elementVertices[5]];
- const int g = __verticesCorrepondance[elementVertices[6]];
- const int h = __verticesCorrepondance[elementVertices[7]];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ const int b = m_mesh_data.__verticesCorrepondance[elementVertices[1]];
+ const int c = m_mesh_data.__verticesCorrepondance[elementVertices[2]];
+ const int d = m_mesh_data.__verticesCorrepondance[elementVertices[3]];
+ const int e = m_mesh_data.__verticesCorrepondance[elementVertices[4]];
+ const int f = m_mesh_data.__verticesCorrepondance[elementVertices[5]];
+ const int g = m_mesh_data.__verticesCorrepondance[elementVertices[6]];
+ const int h = m_mesh_data.__verticesCorrepondance[elementVertices[7]];
if ((a < 0) or (b < 0) or (c < 0) or (d < 0) or (e < 0) or (f < 0) or (g < 0) or (h < 0)) {
throw NormalError("reading file '" + m_filename + "': error reading element " + std::to_string(i) +
" [bad vertices definition]");
}
- __hexahedra[hexahedronNumber] = Hexahedron(a, b, c, d, e, f, g, h);
- __hexahedra_ref[hexahedronNumber] = __references[i];
- __hexahedra_number[hexahedronNumber] = __elementNumber[i];
+ m_mesh_data.__hexahedra[hexahedronNumber] = GmshData::Hexahedron(a, b, c, d, e, f, g, h);
+ m_mesh_data.__hexahedra_ref[hexahedronNumber] = m_mesh_data.__references[i];
+ m_mesh_data.__hexahedra_number[hexahedronNumber] = m_mesh_data.__elementNumber[i];
hexahedronNumber++; // one more hexahedron
break;
}
// Unsupported entities
case 14: { // point
- int& point_number = elementNumber[14];
- const int a = __verticesCorrepondance[elementVertices[0]];
- __points[point_number] = a;
- __points_ref[point_number] = __references[i];
- __points_number[point_number] = __elementNumber[i];
+ int& point_number = elementNumber[14];
+ const int a = m_mesh_data.__verticesCorrepondance[elementVertices[0]];
+ m_mesh_data.__points[point_number] = a;
+ m_mesh_data.__points_ref[point_number] = m_mesh_data.__references[i];
+ m_mesh_data.__points_number[point_number] = m_mesh_data.__elementNumber[i];
point_number++;
break;
}
@@ -629,463 +1200,57 @@ GmshReader::__proceedData()
if ((dimension3_mask, elementNumber) > 0) {
const size_t nb_cells = (dimension3_mask, elementNumber);
- ConnectivityDescriptor descriptor;
-
- descriptor.node_number_vector = __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 = __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] = __tetrahedra[j][r];
- }
- descriptor.cell_type_vector[j] = CellType::Tetrahedron;
- descriptor.cell_number_vector[j] = __tetrahedra_number[j];
- }
- const size_t nb_hexahedra = __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] = __hexahedra[j][r];
- }
- descriptor.cell_type_vector[jh] = CellType::Hexahedron;
- descriptor.cell_number_vector[jh] = __hexahedra_number[j];
- }
-
- std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
- for (unsigned int r = 0; r < __tetrahedra_ref.size(); ++r) {
- const unsigned int elem_number = __tetrahedra_ref[r];
- const unsigned int& ref = __tetrahedra_ref[r];
- ref_cells_map[ref].push_back(elem_number);
- }
-
- for (unsigned int j = 0; j < __hexahedra_ref.size(); ++j) {
- const size_t elem_number = nb_tetrahedra + j;
- const unsigned int& ref = __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_physical_ref_map.at(ref_cell_list.first);
- descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
- }
-
- MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(descriptor);
-
- const auto& node_number_vector = descriptor.node_number_vector;
-
- {
- using Face = ConnectivityFace<3>;
- const std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map = [&] {
- std::unordered_map<Face, FaceId, typename Face::Hash> face_to_id_map;
- for (FaceId l = 0; l < descriptor.face_to_node_vector.size(); ++l) {
- const auto& node_vector = descriptor.face_to_node_vector[l];
- face_to_id_map[Face(node_vector, node_number_vector)] = l;
- }
- return face_to_id_map;
- }();
-
- std::unordered_map<int, FaceId> face_number_id_map = [&] {
- std::unordered_map<int, FaceId> face_number_id_map;
- for (size_t l = 0; l < descriptor.face_number_vector.size(); ++l) {
- face_number_id_map[descriptor.face_number_vector[l]] = l;
- }
- Assert(face_number_id_map.size() == descriptor.face_number_vector.size());
- return face_number_id_map;
- }();
-
- std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
- for (unsigned int f = 0; f < __triangles.size(); ++f) {
- const unsigned int face_id = [&] {
- auto i =
- face_to_id_map.find(Face({__triangles[f][0], __triangles[f][1], __triangles[f][2]}, node_number_vector));
- if (i == face_to_id_map.end()) {
- throw NormalError("face not found");
- }
- return i->second;
- }();
-
- const unsigned int& ref = __triangles_ref[f];
- ref_faces_map[ref].push_back(face_id);
-
- if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
- if (auto i_face = face_number_id_map.find(__quadrangles_number[f]); i_face != face_number_id_map.end()) {
- const int other_face_id = i_face->second;
- std::swap(descriptor.face_number_vector[face_id], descriptor.face_number_vector[other_face_id]);
-
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
-
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
- face_number_id_map[descriptor.face_number_vector[other_face_id]] = other_face_id;
- } else {
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- descriptor.face_number_vector[face_id] = __quadrangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
- }
- }
- }
-
- for (unsigned int f = 0; f < __quadrangles.size(); ++f) {
- const unsigned int face_id = [&] {
- auto i = face_to_id_map.find(
- Face({__quadrangles[f][0], __quadrangles[f][1], __quadrangles[f][2], __quadrangles[f][3]},
- node_number_vector));
- if (i == face_to_id_map.end()) {
- throw NormalError("face not found");
- }
- return i->second;
- }();
-
- const unsigned int& ref = __quadrangles_ref[f];
- ref_faces_map[ref].push_back(face_id);
-
- if (descriptor.face_number_vector[face_id] != __quadrangles_number[f]) {
- if (auto i_face = face_number_id_map.find(__quadrangles_number[f]); i_face != face_number_id_map.end()) {
- const int other_face_id = i_face->second;
- std::swap(descriptor.face_number_vector[face_id], descriptor.face_number_vector[other_face_id]);
-
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- face_number_id_map.erase(descriptor.face_number_vector[other_face_id]);
-
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
- face_number_id_map[descriptor.face_number_vector[other_face_id]] = other_face_id;
- } else {
- face_number_id_map.erase(descriptor.face_number_vector[face_id]);
- descriptor.face_number_vector[face_id] = __quadrangles_number[f];
- face_number_id_map[descriptor.face_number_vector[face_id]] = face_id;
- }
- }
- }
-
- for (const auto& ref_face_list : ref_faces_map) {
- Array<FaceId> face_list(ref_face_list.second.size());
- for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
- face_list[j] = ref_face_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
- descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
- }
- }
- MeshBuilderBase::_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 < __edges.size(); ++e) {
- const unsigned int edge_id = [&] {
- auto i = edge_to_id_map.find(Edge({__edges[e][0], __edges[e][1]}, node_number_vector));
- if (i == edge_to_id_map.end()) {
- std::stringstream error_msg;
- error_msg << "edge " << __edges[e][0] << " not found";
- throw NormalError(error_msg.str());
- }
- return i->second;
- }();
- const unsigned int& ref = __edges_ref[e];
- ref_edges_map[ref].push_back(edge_id);
-
- if (descriptor.edge_number_vector[edge_id] != __edges_number[e]) {
- if (auto i_edge = edge_number_id_map.find(__edges_number[e]); i_edge != edge_number_id_map.end()) {
- const int other_edge_id = i_edge->second;
- std::swap(descriptor.edge_number_vector[edge_id], descriptor.edge_number_vector[other_edge_id]);
-
- edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
- edge_number_id_map.erase(descriptor.edge_number_vector[other_edge_id]);
-
- edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_id;
- edge_number_id_map[descriptor.edge_number_vector[other_edge_id]] = other_edge_id;
- } else {
- edge_number_id_map.erase(descriptor.edge_number_vector[edge_id]);
- descriptor.edge_number_vector[edge_id] = __edges_number[e];
- edge_number_id_map[descriptor.edge_number_vector[edge_id]] = edge_id;
- }
- }
- }
-
- for (const auto& ref_edge_list : ref_edges_map) {
- Array<EdgeId> edge_list(ref_edge_list.second.size());
- for (size_t j = 0; j < ref_edge_list.second.size(); ++j) {
- edge_list[j] = ref_edge_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_edge_list.first);
- descriptor.addRefItemList(RefEdgeList{physical_ref_id.refId(), edge_list});
- }
- }
-
- std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
- for (unsigned int r = 0; r < __points.size(); ++r) {
- const unsigned int point_number = __points[r];
- const unsigned int& ref = __points_ref[r];
- ref_points_map[ref].push_back(point_number);
- }
-
- for (const auto& ref_point_list : ref_points_map) {
- Array<NodeId> point_list(ref_point_list.second.size());
- for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
- point_list[j] = ref_point_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_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());
+ GmshConnectivityBuilder<3> connectivity_builder(m_mesh_data, nb_cells);
- std::shared_ptr p_connectivity = Connectivity3D::build(descriptor);
- Connectivity3D& connectivity = *p_connectivity;
+ std::shared_ptr p_connectivity =
+ std::dynamic_pointer_cast<const Connectivity3D>(connectivity_builder.connectivity());
+ const Connectivity3D& connectivity = *p_connectivity;
using MeshType = Mesh<Connectivity3D>;
using Rd = TinyVector<3, double>;
NodeValue<Rd> xr(connectivity);
- for (NodeId i = 0; i < __vertices.size(); ++i) {
- xr[i][0] = __vertices[i][0];
- xr[i][1] = __vertices[i][1];
- xr[i][2] = __vertices[i][2];
+ for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+ xr[i][0] = m_mesh_data.__vertices[i][0];
+ xr[i][1] = m_mesh_data.__vertices[i][1];
+ xr[i][2] = m_mesh_data.__vertices[i][2];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
} else if ((dimension2_mask, elementNumber) > 0) {
const size_t nb_cells = (dimension2_mask, elementNumber);
- ConnectivityDescriptor descriptor;
-
- descriptor.node_number_vector = __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_triangles = __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] = __triangles[j][r];
- }
- descriptor.cell_type_vector[j] = CellType::Triangle;
- descriptor.cell_number_vector[j] = __triangles_number[j];
- }
-
- const size_t nb_quadrangles = __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] = __quadrangles[j][r];
- }
- descriptor.cell_type_vector[jq] = CellType::Quadrangle;
- descriptor.cell_number_vector[jq] = __quadrangles_number[j];
- }
-
- std::map<unsigned int, std::vector<unsigned int>> ref_cells_map;
- for (unsigned int r = 0; r < __triangles_ref.size(); ++r) {
- const unsigned int elem_number = __triangles_ref[r];
- const unsigned int& ref = __triangles_ref[r];
- ref_cells_map[ref].push_back(elem_number);
- }
-
- for (unsigned int j = 0; j < __quadrangles_ref.size(); ++j) {
- const size_t elem_number = nb_triangles + j;
- const unsigned int& ref = __quadrangles_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_physical_ref_map.at(ref_cell_list.first);
- descriptor.addRefItemList(RefCellList(physical_ref_id.refId(), cell_list));
- }
-
- MeshBuilderBase::_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 = [&] {
- 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 e = 0; e < __edges.size(); ++e) {
- const unsigned int edge_id = [&] {
- auto i = face_to_id_map.find(Face({__edges[e][0], __edges[e][1]}, node_number_vector));
- if (i == face_to_id_map.end()) {
- std::stringstream error_msg;
- error_msg << "face " << __edges[e][0] << " not found";
- throw NormalError(error_msg.str());
- }
- return i->second;
- }();
- const unsigned int& ref = __edges_ref[e];
- ref_faces_map[ref].push_back(edge_id);
-
- if (descriptor.face_number_vector[edge_id] != __edges_number[e]) {
- if (auto i_face = face_number_id_map.find(__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]);
-
- 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[descriptor.face_number_vector[edge_id]] = edge_id;
- face_number_id_map[descriptor.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] = __edges_number[e];
- face_number_id_map[descriptor.face_number_vector[edge_id]] = edge_id;
- }
- }
- }
-
- for (const auto& ref_face_list : ref_faces_map) {
- Array<FaceId> face_list(ref_face_list.second.size());
- for (size_t j = 0; j < ref_face_list.second.size(); ++j) {
- face_list[j] = ref_face_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_face_list.first);
- descriptor.addRefItemList(RefFaceList{physical_ref_id.refId(), face_list});
- }
-
- std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
- for (unsigned int r = 0; r < __points.size(); ++r) {
- const unsigned int point_number = __points[r];
- const unsigned int& ref = __points_ref[r];
- ref_points_map[ref].push_back(point_number);
- }
-
- for (const auto& ref_point_list : ref_points_map) {
- Array<NodeId> point_list(ref_point_list.second.size());
- for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
- point_list[j] = ref_point_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_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.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ GmshConnectivityBuilder<2> connectivity_builder(m_mesh_data, nb_cells);
- std::shared_ptr p_connectivity = Connectivity2D::build(descriptor);
- Connectivity2D& connectivity = *p_connectivity;
+ std::shared_ptr p_connectivity =
+ std::dynamic_pointer_cast<const Connectivity2D>(connectivity_builder.connectivity());
+ const Connectivity2D& connectivity = *p_connectivity;
using MeshType = Mesh<Connectivity2D>;
using Rd = TinyVector<2, double>;
NodeValue<Rd> xr(connectivity);
- for (NodeId i = 0; i < __vertices.size(); ++i) {
- xr[i][0] = __vertices[i][0];
- xr[i][1] = __vertices[i][1];
+ for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+ xr[i][0] = m_mesh_data.__vertices[i][0];
+ xr[i][1] = m_mesh_data.__vertices[i][1];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
} else if ((dimension1_mask, elementNumber) > 0) {
const size_t nb_cells = (dimension1_mask, elementNumber);
- ConnectivityDescriptor descriptor;
-
- descriptor.node_number_vector = __verticesNumbers;
- descriptor.cell_type_vector.resize(nb_cells);
- descriptor.cell_number_vector.resize(nb_cells);
- descriptor.cell_to_node_vector.resize(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] = __edges[j][r];
- }
- descriptor.cell_type_vector[j] = CellType::Line;
- descriptor.cell_number_vector[j] = __edges_number[j];
- }
-
- std::map<unsigned int, std::vector<unsigned int>> ref_points_map;
- for (unsigned int r = 0; r < __points.size(); ++r) {
- const unsigned int point_number = __points[r];
- const unsigned int& ref = __points_ref[r];
- ref_points_map[ref].push_back(point_number);
- }
-
- for (const auto& ref_point_list : ref_points_map) {
- Array<NodeId> point_list(ref_point_list.second.size());
- for (size_t j = 0; j < ref_point_list.second.size(); ++j) {
- point_list[j] = ref_point_list.second[j];
- }
- const PhysicalRefId& physical_ref_id = m_physical_ref_map.at(ref_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.node_owner_vector.resize(descriptor.node_number_vector.size());
- std::fill(descriptor.node_owner_vector.begin(), descriptor.node_owner_vector.end(), parallel::rank());
+ GmshConnectivityBuilder<1> connectivity_builder(m_mesh_data, nb_cells);
- std::shared_ptr p_connectivity = Connectivity1D::build(descriptor);
- Connectivity1D& connectivity = *p_connectivity;
+ std::shared_ptr p_connectivity =
+ std::dynamic_pointer_cast<const Connectivity1D>(connectivity_builder.connectivity());
+ const Connectivity1D& connectivity = *p_connectivity;
using MeshType = Mesh<Connectivity1D>;
using Rd = TinyVector<1, double>;
NodeValue<Rd> xr(connectivity);
- for (NodeId i = 0; i < __vertices.size(); ++i) {
- xr[i][0] = __vertices[i][0];
+ for (NodeId i = 0; i < m_mesh_data.__vertices.size(); ++i) {
+ xr[i][0] = m_mesh_data.__vertices[i][0];
}
m_mesh = std::make_shared<MeshType>(p_connectivity, xr);
diff --git a/src/mesh/GmshReader.hpp b/src/mesh/GmshReader.hpp
index 4edd1fc7b..ae18d30ca 100644
--- a/src/mesh/GmshReader.hpp
+++ b/src/mesh/GmshReader.hpp
@@ -69,73 +69,80 @@ class GmshReader : public MeshBuilderBase
~PhysicalRefId() = default;
};
+ struct GmshData
+ {
+ /**
+ * Gmsh format provides a numbered, none ordrered and none dense
+ * vertices list, this stores the number of read vertices.
+ */
+ std::vector<int> __verticesNumbers;
+
+ Array<R3> __vertices;
+
+ using Point = unsigned int;
+ Array<Point> __points;
+ std::vector<int> __points_ref;
+ std::vector<int> __points_number;
+
+ using Edge = TinyVector<2, unsigned int>;
+ Array<Edge> __edges;
+ std::vector<int> __edges_ref;
+ std::vector<int> __edges_number;
+
+ using Triangle = TinyVector<3, unsigned int>;
+ Array<Triangle> __triangles;
+ std::vector<int> __triangles_ref;
+ std::vector<int> __triangles_number;
+
+ using Quadrangle = TinyVector<4, unsigned int>;
+ Array<Quadrangle> __quadrangles;
+ std::vector<int> __quadrangles_ref;
+ std::vector<int> __quadrangles_number;
+
+ using Tetrahedron = TinyVector<4, unsigned int>;
+ Array<Tetrahedron> __tetrahedra;
+ std::vector<int> __tetrahedra_ref;
+ std::vector<int> __tetrahedra_number;
+
+ using Hexahedron = TinyVector<8, unsigned int>;
+ Array<Hexahedron> __hexahedra;
+ std::vector<int> __hexahedra_ref;
+ std::vector<int> __hexahedra_number;
+
+ /**
+ * Gmsh format provides a numbered, none ordrered and none dense
+ * vertices list, this provides vertices renumbering correspondence
+ */
+ std::vector<int> __verticesCorrepondance;
+
+ /**
+ * elements types
+ */
+ std::vector<int> __elementNumber;
+
+ /**
+ * elements types
+ */
+ std::vector<short> __elementType;
+
+ /**
+ * References
+ */
+ std::vector<int> __references;
+
+ /**
+ * References
+ */
+ std::vector<std::vector<int> > __elementVertices;
+
+ std::map<unsigned int, PhysicalRefId> m_physical_ref_map;
+ };
+
private:
std::ifstream m_fin;
const std::string m_filename;
- /**
- * Gmsh format provides a numbered, none ordrered and none dense
- * vertices list, this stores the number of read vertices.
- */
- std::vector<int> __verticesNumbers;
-
- Array<R3> __vertices;
-
- using Point = unsigned int;
- Array<Point> __points;
- std::vector<int> __points_ref;
- std::vector<int> __points_number;
-
- using Edge = TinyVector<2, unsigned int>;
- Array<Edge> __edges;
- std::vector<int> __edges_ref;
- std::vector<int> __edges_number;
-
- using Triangle = TinyVector<3, unsigned int>;
- Array<Triangle> __triangles;
- std::vector<int> __triangles_ref;
- std::vector<int> __triangles_number;
-
- using Quadrangle = TinyVector<4, unsigned int>;
- Array<Quadrangle> __quadrangles;
- std::vector<int> __quadrangles_ref;
- std::vector<int> __quadrangles_number;
-
- using Tetrahedron = TinyVector<4, unsigned int>;
- Array<Tetrahedron> __tetrahedra;
- std::vector<int> __tetrahedra_ref;
- std::vector<int> __tetrahedra_number;
-
- using Hexahedron = TinyVector<8, unsigned int>;
- Array<Hexahedron> __hexahedra;
- std::vector<int> __hexahedra_ref;
- std::vector<int> __hexahedra_number;
-
- /**
- * Gmsh format provides a numbered, none ordrered and none dense
- * vertices list, this provides vertices renumbering correspondance
- */
- std::vector<int> __verticesCorrepondance;
-
- /**
- * elements types
- */
- std::vector<int> __elementNumber;
-
- /**
- * elements types
- */
- std::vector<short> __elementType;
-
- /**
- * References
- */
- std::vector<int> __references;
-
- /**
- * References
- */
- std::vector<std::vector<int> > __elementVertices;
+ GmshData m_mesh_data;
/**
* Stores the number of nodes associated to each primitive
@@ -155,11 +162,6 @@ class GmshReader : public MeshBuilderBase
*/
std::map<int, std::string> __primitivesNames;
- bool __binary; /**< true if binary format */
- bool __convertEndian; /**< true if needs to adapt endianess */
-
- std::map<unsigned int, PhysicalRefId> m_physical_ref_map;
-
int
_getInteger()
{
diff --git a/src/mesh/MeshBuilderBase.cpp b/src/mesh/MeshBuilderBase.cpp
index 7f7221b04..fc9020325 100644
--- a/src/mesh/MeshBuilderBase.cpp
+++ b/src/mesh/MeshBuilderBase.cpp
@@ -38,9 +38,15 @@ MeshBuilderBase::_dispatch()
m_mesh = std::make_shared<MeshType>(dispatched_connectivity, dispatched_xr);
}
+template void MeshBuilderBase::_dispatch<1>();
+template void MeshBuilderBase::_dispatch<2>();
+template void MeshBuilderBase::_dispatch<3>();
+
+/// ============
+
template <size_t Dimension>
void
-MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor)
+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>;
@@ -247,7 +253,7 @@ MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities(ConnectivityDescripto
template <size_t Dimension>
void
-MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
+ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor)
{
static_assert(Dimension == 3, "Invalid dimension to compute face-edge connectivities");
using FaceEdgeInfo = std::tuple<FaceId, unsigned short, bool>;
@@ -437,12 +443,8 @@ MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(Connectivi
}
}
-template void MeshBuilderBase::_dispatch<1>();
-template void MeshBuilderBase::_dispatch<2>();
-template void MeshBuilderBase::_dispatch<3>();
-
-template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
-template void MeshBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<2>(ConnectivityDescriptor& descriptor);
+template void ConnectivityBuilderBase::_computeCellFaceAndFaceNodeConnectivities<3>(ConnectivityDescriptor& descriptor);
-template void MeshBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
+template void ConnectivityBuilderBase::_computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities<3>(
ConnectivityDescriptor& descriptor);
diff --git a/src/mesh/MeshBuilderBase.hpp b/src/mesh/MeshBuilderBase.hpp
index 03eae11d5..24bd05582 100644
--- a/src/mesh/MeshBuilderBase.hpp
+++ b/src/mesh/MeshBuilderBase.hpp
@@ -2,6 +2,9 @@
#define MESH_BUILDER_BASE_HPP
#include <mesh/IMesh.hpp>
+
+#include <mesh/IConnectivity.hpp>
+
#include <mesh/ItemId.hpp>
#include <utils/PugsAssert.hpp>
#include <utils/PugsMacros.hpp>
@@ -11,13 +14,13 @@
class ConnectivityDescriptor;
-class MeshBuilderBase
+class ConnectivityBuilderBase
{
protected:
template <size_t Dimension>
class ConnectivityFace;
- std::shared_ptr<const IMesh> m_mesh;
+ std::shared_ptr<const IConnectivity> m_connectivity;
template <size_t Dimension>
static void _computeCellFaceAndFaceNodeConnectivities(ConnectivityDescriptor& descriptor);
@@ -25,22 +28,19 @@ class MeshBuilderBase
template <size_t Dimension>
static void _computeFaceEdgeAndEdgeNodeAndCellEdgeConnectivities(ConnectivityDescriptor& descriptor);
- template <int Dimension>
- void _dispatch();
-
public:
- std::shared_ptr<const IMesh>
- mesh() const
+ std::shared_ptr<const IConnectivity>
+ connectivity() const
{
- return m_mesh;
+ return m_connectivity;
}
- MeshBuilderBase() = default;
- ~MeshBuilderBase() = default;
+ ConnectivityBuilderBase() = default;
+ ~ConnectivityBuilderBase() = default;
};
template <>
-class MeshBuilderBase::ConnectivityFace<2>
+class ConnectivityBuilderBase::ConnectivityFace<2>
{
public:
friend struct Hash;
@@ -119,7 +119,7 @@ class MeshBuilderBase::ConnectivityFace<2>
};
template <>
-class MeshBuilderBase::ConnectivityFace<3>
+class ConnectivityBuilderBase::ConnectivityFace<3>
{
public:
friend struct Hash;
@@ -225,4 +225,23 @@ class MeshBuilderBase::ConnectivityFace<3>
~ConnectivityFace() = default;
};
+class MeshBuilderBase
+{
+ protected:
+ std::shared_ptr<const IMesh> m_mesh;
+
+ template <int Dimension>
+ void _dispatch();
+
+ public:
+ std::shared_ptr<const IMesh>
+ mesh() const
+ {
+ return m_mesh;
+ }
+
+ MeshBuilderBase() = default;
+ ~MeshBuilderBase() = default;
+};
+
#endif // MESH_BUILDER_BASE_HPP
--
GitLab