diff --git a/src/main.cpp b/src/main.cpp
index a0ab842a3ca1b0961ce64fddfeb1fb04d68cf714..e351d6471a50bb133286a1b53bbd9bc7ee12f1dc 100644
--- a/src/main.cpp
+++ b/src/main.cpp
@@ -122,7 +122,9 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t);
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}},t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -136,7 +138,9 @@ int main(int argc, char *argv[])
}
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
@@ -233,7 +237,9 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t);
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}}, t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -247,7 +253,9 @@ int main(int argc, char *argv[])
}
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
@@ -333,7 +341,9 @@ int main(int argc, char *argv[])
while((t<tmax) and (iteration<itermax)) {
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t);
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}}, t);
double dt = 0.4*acoustic_solver.acoustic_dt(Vj, cj);
if (t+dt>tmax) {
dt=tmax-t;
@@ -346,7 +356,9 @@ int main(int argc, char *argv[])
}
vtk_writer.write(mesh, {NamedItemValue{"density", rhoj},
NamedItemValue{"velocity", unknowns.uj()},
- NamedItemValue{"coords", mesh.xr()}}, t, true); // forces last output
+ NamedItemValue{"coords", mesh.xr()},
+ NamedItemValue{"cell_owner", mesh.connectivity().cellOwner()},
+ NamedItemValue{"node_owner", mesh.connectivity().nodeOwner()}}, t, true); // forces last output
pout() << "* " << rang::style::underline << "Final time" << rang::style::reset
<< ": " << rang::fgB::green << t << rang::fg::reset << " (" << iteration << " iterations)\n";
diff --git a/src/mesh/Connectivity.cpp b/src/mesh/Connectivity.cpp
index bc8e7093aabf110b3bcebab71fddd5f2585883e5..07701bbb950f6d1cf6ecccff9982b54a86146b67 100644
--- a/src/mesh/Connectivity.cpp
+++ b/src/mesh/Connectivity.cpp
@@ -218,7 +218,9 @@ Connectivity<Dimension>::
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector,
- const std::vector<int>& node_number_vector)
+ const std::vector<int>& node_number_vector,
+ const std::vector<int>& cell_owner_vector,
+ const std::vector<int>& node_owner_vector)
{
Assert(cell_by_node_vector.size() == cell_type_vector.size());
Assert(cell_number_vector.size() == cell_type_vector.size());
@@ -268,6 +270,18 @@ Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
m_inv_cell_nb_nodes = inv_cell_nb_nodes;
}
+ {
+ CellValue<int> cell_owner(*this);
+ cell_owner = convert_to_array(cell_owner_vector);
+ m_cell_owner = cell_owner;
+ }
+
+ {
+ NodeValue<int> node_owner(*this);
+ node_owner = convert_to_array(node_owner_vector);
+ m_node_owner = node_owner;
+ }
+
if constexpr (Dimension>1) {
this->_computeCellFaceAndFaceNodeConnectivities();
}
@@ -278,16 +292,22 @@ template Connectivity1D::
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector,
- const std::vector<int>& node_number_vector);
+ const std::vector<int>& node_number_vector,
+ const std::vector<int>& cell_owner_vector,
+ const std::vector<int>& node_owner_vector);
template Connectivity2D::
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector,
- const std::vector<int>& node_number_vector);
+ const std::vector<int>& node_number_vector,
+ const std::vector<int>& cell_owner_vector,
+ const std::vector<int>& node_owner_vector);
template Connectivity3D::
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector,
- const std::vector<int>& node_number_vector);
+ const std::vector<int>& node_number_vector,
+ const std::vector<int>& cell_owner_vector,
+ const std::vector<int>& node_owner_vector);
diff --git a/src/mesh/Connectivity.hpp b/src/mesh/Connectivity.hpp
index c786eb6c316ef8f110e2f5521378fc829f5ce35b..b8783f321c82f0eaa8d5c215b4dec0aefbe9c460 100644
--- a/src/mesh/Connectivity.hpp
+++ b/src/mesh/Connectivity.hpp
@@ -251,6 +251,9 @@ class Connectivity final
NodeValue<const int> m_node_number;
+ CellValue<const int> m_cell_owner;
+ NodeValue<const int> m_node_owner;
+
FaceValuePerCell<const bool> m_cell_face_is_reversed;
NodeValuePerCell<const unsigned short> m_cell_local_numbers_in_their_nodes;
@@ -332,6 +335,18 @@ class Connectivity final
return m_node_number;
}
+ PASTIS_INLINE
+ const CellValue<const int>& cellOwner() const
+ {
+ return m_cell_owner;
+ }
+
+ PASTIS_INLINE
+ const NodeValue<const int>& nodeOwner() const
+ {
+ return m_node_owner;
+ }
+
PASTIS_INLINE
const bool& isConnectivityMatrixBuilt(const ItemType& item_type_0,
const ItemType& item_type_1) const
@@ -662,7 +677,9 @@ class Connectivity final
Connectivity(const std::vector<std::vector<unsigned int>>& cell_by_node_vector,
const std::vector<CellType>& cell_type_vector,
const std::vector<int>& cell_number_vector,
- const std::vector<int>& node_number_vector);
+ const std::vector<int>& node_number_vector,
+ const std::vector<int>& cell_owner_vector,
+ const std::vector<int>& node_owner_vector);
~Connectivity()
{
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 4ec88d09c1b0dfe091be78282e9784e7885792d8..6a19ba4f2595e5b7c7bb9f5f6ae78e9c9607ce07 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -33,11 +33,11 @@ class ErrorHandler
};
private:
- const std::string __filename; /**< The source file name where the error occured */
- const size_t __lineNumber; /**< The line number where exception was raised */
+ const std::string __filename; /**< The source file name where the error occured */
+ const size_t __lineNumber; /**< The line number where exception was raised */
const std::string __errorMessage; /**< The reporting message */
- const Type __type; /**< the type of the error */
+ const Type __type; /**< the type of the error */
public:
/**
* Prints the error message
@@ -135,23 +135,56 @@ class MeshDispatcher
private:
const MeshType& m_mesh;
+ CellValue<const int> m_cell_new_owner;
+ NodeValue<const int> m_node_new_owner;
+
using CellListToSendByProc = std::vector<Array<const CellId>>;
const CellListToSendByProc m_cell_list_to_send_by_proc;
Array<int> m_nb_cell_to_send_by_proc;
Array<int> m_nb_cell_to_recv_by_proc;
- const CellListToSendByProc
- _buildCellListToSend() const
+ CellValue<int> _getCellNewOwner()
{
CSRGraph mesh_graph = m_mesh.cellToCellGraph();
-
Partitioner P;
- Array<int> cell_new_owner = P.partition(mesh_graph);
+ CellValue<int> cell_new_owner(m_mesh.connectivity());
+ cell_new_owner = P.partition(mesh_graph);
+ return cell_new_owner;
+ }
+
+ NodeValue<int> _getNodeNewOwner()
+ {
+ const auto& node_to_cell_matrix
+ = m_mesh.connectivity().nodeToCellMatrix();
+
+ pout() << __FILE__ << ':' << __LINE__ << ": use Min function\n";
+#warning could use a better policy
+ NodeValue<int> node_new_owner(m_mesh.connectivity());
+ parallel_for(m_mesh.numberOfNodes(), PASTIS_LAMBDA(const NodeId& r) {
+ const auto& node_to_cell = node_to_cell_matrix[r];
+ CellId Jmin = node_to_cell[0];
+
+ for (size_t j=1; j<node_to_cell.size(); ++j) {
+ const CellId J = node_to_cell[j];
+ if (J<Jmin) {
+ Jmin=J;
+ }
+ }
+ node_new_owner[r] = m_cell_new_owner[Jmin];
+ });
+
+#warning Add missing synchronize
+ return node_new_owner;
+ }
+
+ const CellListToSendByProc
+ _buildCellListToSend() const
+ {
std::vector<std::vector<CellId>> cell_vector_to_send_by_proc(parallel::size());
- for (size_t i=0; i<cell_new_owner.size(); ++i) {
- cell_vector_to_send_by_proc[cell_new_owner[i]].push_back(CellId{i});
+ for (CellId j=0; j<m_mesh.numberOfCells(); ++j) {
+ cell_vector_to_send_by_proc[m_cell_new_owner[j]].push_back(j);
}
std::vector<Array<const CellId>> cell_list_to_send_by_proc(parallel::size());
@@ -173,6 +206,18 @@ class MeshDispatcher
public:
+ PASTIS_INLINE
+ const CellValue<const int>& cellNewOwner() const
+ {
+ return m_cell_new_owner;
+ }
+
+ PASTIS_INLINE
+ const NodeValue<const int>& nodeNewOwner() const
+ {
+ return m_node_new_owner;
+ }
+
template<typename DataType>
std::vector<Array<std::remove_const_t<DataType>>>
exchange(const CellValue<DataType>& cell_value) const
@@ -205,6 +250,8 @@ class MeshDispatcher
MeshDispatcher(const MeshType& mesh)
: m_mesh(mesh),
+ m_cell_new_owner(_getCellNewOwner()),
+ m_node_new_owner(_getNodeNewOwner()),
m_cell_list_to_send_by_proc(_buildCellListToSend()),
m_nb_cell_to_send_by_proc(_buildNbCellToSend()),
m_nb_cell_to_recv_by_proc(parallel::allToAll(m_nb_cell_to_send_by_proc))
@@ -226,7 +273,7 @@ void GmshReader::_dispatch()
using MeshType = Mesh<ConnectivityType>;
if (not m_mesh) {
- std::shared_ptr<ConnectivityType> connectivity(new ConnectivityType({},{},{},{}));
+ std::shared_ptr<ConnectivityType> connectivity(new ConnectivityType({},{},{},{}, {}, {}));
NodeValue<Rd> xr;
m_mesh = std::shared_ptr<IMesh>(new MeshType(connectivity, xr));
}
@@ -234,13 +281,15 @@ void GmshReader::_dispatch()
MeshDispatcher<Dimension> dispatcher(mesh);
-
std::vector<Array<int>> recv_cell_number_by_proc
= dispatcher.exchange(mesh.connectivity().cellNumber());
std::vector<Array<CellType>> recv_cell_type_by_proc
= dispatcher.exchange(mesh.connectivity().cellType());
+ std::vector<Array<int>> recv_cell_new_owner_by_proc
+ = dispatcher.exchange(dispatcher.cellNewOwner());
+
CellValue<int> number_of_node_per_cell(mesh.connectivity());
const auto& cell_to_node_matrix = mesh.connectivity().cellToNodeMatrix();
@@ -477,11 +526,50 @@ void GmshReader::_dispatch()
using ConnectivityType = Connectivity<Dimension>;
+ std::vector<int> cell_new_owner_vector(cell_number_id_map.size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ for (size_t i=0; i<recv_cell_number_by_proc[i_rank].size(); ++i) {
+ int cell_number = recv_cell_number_by_proc[i_rank][i];
+ const auto& searched_cell_id = cell_number_id_map.find(cell_number);
+ Assert(searched_cell_id != cell_number_id_map.end());
+ cell_new_owner_vector[searched_cell_id->second] = recv_cell_new_owner_by_proc[i_rank][i];
+ }
+ }
+
+
+ const NodeValue<const int>& new_node_owner = dispatcher.nodeNewOwner();
+ std::vector<Array<const int>> send_node_owner_by_proc(parallel::size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ Array<int> send_node_owner(nb_node_to_send_by_proc[i_rank]);
+ const Array<const NodeId> send_node_id = send_node_id_by_proc[i_rank];
+ parallel_for(send_node_id.size(), PASTIS_LAMBDA(const size_t& r) {
+ const NodeId& node_id = send_node_id[r];
+ send_node_owner[r] = new_node_owner[node_id];
+ });
+ send_node_owner_by_proc[i_rank] = send_node_owner;
+ }
+
+ std::vector<Array<int>> recv_node_owner_by_proc(parallel::size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ recv_node_owner_by_proc[i_rank] = Array<int>(nb_node_to_recv_by_proc[i_rank]);
+ }
+ parallel::exchange(send_node_owner_by_proc, recv_node_owner_by_proc);
+
+ std::vector<int> node_new_owner_vector(new_node_number.size());
+ for (size_t i_rank=0; i_rank<parallel::size(); ++i_rank) {
+ for (size_t r=0; r<recv_node_owner_by_proc[i_rank].size(); ++r) {
+ const NodeId& node_id = recv_node_id_correspondance_by_proc[i_rank][r];
+ node_new_owner_vector[node_id] = recv_node_owner_by_proc[i_rank][r];
+ }
+ }
+
std::shared_ptr<ConnectivityType> p_connectivity
= std::make_shared<ConnectivityType>(cell_by_node_vector,
new_cell_type,
new_cell_number,
- new_node_number);
+ new_node_number,
+ cell_new_owner_vector,
+ node_new_owner_vector);
const NodeValue<const Rd>& xr = mesh.xr();
std::vector<Array<const Rd>> send_node_coord_by_proc(parallel::size());
@@ -1253,10 +1341,18 @@ GmshReader::__proceedData()
cell_number_vector[jh] = __hexahedra_number[j];
}
+ std::vector<int> cell_owner_vector(nb_cells);
+ std::fill(cell_owner_vector.begin(), cell_owner_vector.end(), parallel::rank());
+ std::vector<int> node_owner_vector(node_number_vector.size());
+ std::fill(node_owner_vector.begin(), node_owner_vector.end(), parallel::rank());
+
+
std::shared_ptr<Connectivity3D> p_connectivity(new Connectivity3D(cell_by_node_vector,
cell_type_vector,
cell_number_vector,
- node_number_vector));
+ node_number_vector,
+ cell_owner_vector,
+ node_owner_vector));
Connectivity3D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -1321,10 +1417,17 @@ GmshReader::__proceedData()
cell_number_vector[jq] = __quadrangles_number[j];
}
+ std::vector<int> cell_owner_vector(nb_cells);
+ std::fill(cell_owner_vector.begin(), cell_owner_vector.end(), parallel::rank());
+ std::vector<int> node_owner_vector(node_number_vector.size());
+ std::fill(node_owner_vector.begin(), node_owner_vector.end(), parallel::rank());
+
std::shared_ptr<Connectivity2D> p_connectivity(new Connectivity2D(cell_by_node_vector,
cell_type_vector,
cell_number_vector,
- node_number_vector));
+ node_number_vector,
+ cell_owner_vector,
+ node_owner_vector));
Connectivity2D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_faces_map;
@@ -1385,10 +1488,17 @@ GmshReader::__proceedData()
cell_number_vector[j] = __edges_number[j];
}
+ std::vector<int> cell_owner_vector(nb_cells);
+ std::fill(cell_owner_vector.begin(), cell_owner_vector.end(), parallel::rank());
+ std::vector<int> node_owner_vector(node_number_vector.size());
+ std::fill(node_owner_vector.begin(), node_owner_vector.end(), parallel::rank());
+
std::shared_ptr<Connectivity1D> p_connectivity(new Connectivity1D(cell_by_node_vector,
cell_type_vector,
cell_number_vector,
- node_number_vector));
+ node_number_vector,
+ cell_owner_vector,
+ node_owner_vector));
Connectivity1D& connectivity = *p_connectivity;
std::map<unsigned int, std::vector<unsigned int>> ref_points_map;