diff --git a/src/mesh/CMakeLists.txt b/src/mesh/CMakeLists.txt
index c28b0643fd8bbf1b2677512423d2185d9a7ff970..10ea744770b20545201c5eaa9b6b947de7bd3ee8 100644
--- a/src/mesh/CMakeLists.txt
+++ b/src/mesh/CMakeLists.txt
@@ -9,7 +9,7 @@ add_library(
ConnectivityComputer.cpp
GmshReader.cpp)
-#include_directories(${PASTIS_SOURCE_DIR}/utils)
+include_directories(${PASTIS_SOURCE_DIR}/utils)
# Additional dependencies
#add_dependencies(PastisMesh)
diff --git a/src/mesh/GmshReader.cpp b/src/mesh/GmshReader.cpp
index 3a2d5c6ecc290a6ae71c234dd0551c7b48a22dca..e599ea8083c4d0a25471fab73c22192c1478af50 100644
--- a/src/mesh/GmshReader.cpp
+++ b/src/mesh/GmshReader.cpp
@@ -12,6 +12,8 @@
#include <Mesh.hpp>
#include <RefFaceList.hpp>
+#include <Messenger.hpp>
+#include <Partitioner.hpp>
#include <map>
#include <regex>
@@ -141,162 +143,178 @@ ErrorHandler(const std::string& filename,
GmshReader::GmshReader(const std::string& filename)
: m_filename(filename)
{
- try {
- m_fin.open(m_filename);
- if (not m_fin) {
- perr() << "cannot read file '" << m_filename << "'\n";
- std::exit(0);
- }
+ if (commRank() == 0) {
+ try {
+ m_fin.open(m_filename);
+ if (not m_fin) {
+ perr() << "cannot read file '" << m_filename << "'\n";
+ std::exit(0);
+ }
- // gmsh 2.2 format keywords
- __keywordList["$MeshFormat"] = MESHFORMAT;
- __keywordList["$EndMeshFormat"] = ENDMESHFORMAT;
- __keywordList["$Nodes"] = NODES;
- __keywordList["$EndNodes"] = ENDNODES;
- __keywordList["$Elements"] = ELEMENTS;
- __keywordList["$EndElements"] = ENDELEMENTS;
- __keywordList["$PhysicalNames"] = PHYSICALNAMES;
- __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
-
- __numberOfPrimitiveNodes.resize(16);
- __numberOfPrimitiveNodes[ 0] = 2; // edge
- __numberOfPrimitiveNodes[ 1] = 3; // triangle
- __numberOfPrimitiveNodes[ 2] = 4; // quadrangle
- __numberOfPrimitiveNodes[ 3] = 4; // Tetrahedron
- __numberOfPrimitiveNodes[ 4] = 8; // Hexaredron
- __numberOfPrimitiveNodes[ 5] = 6; // Prism
- __numberOfPrimitiveNodes[ 6] = 5; // Pyramid
- __numberOfPrimitiveNodes[ 7] = 3; // second order edge
- __numberOfPrimitiveNodes[ 8] = 6; // second order triangle
- __numberOfPrimitiveNodes[ 9] = 9; // second order quadrangle
- __numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
- __numberOfPrimitiveNodes[11] = 27; // second order hexahedron
- __numberOfPrimitiveNodes[12] = 18; // second order prism
- __numberOfPrimitiveNodes[13] = 14; // second order pyramid
- __numberOfPrimitiveNodes[14] = 1; // point
-
- __primitivesNames[0] = "edges";
- __supportedPrimitives[0] = true;
- __primitivesNames[1] = "triangles";
- __supportedPrimitives[1] = true;
- __primitivesNames[2] = "quadrangles";
- __supportedPrimitives[2] = true;
- __primitivesNames[3] = "tetrahedra";
- __supportedPrimitives[3] = true;
- __primitivesNames[4] = "hexahedra";
- __supportedPrimitives[4] = true;
- __primitivesNames[5] = "prisms";
- __supportedPrimitives[5] = false;
- __primitivesNames[6] = "pyramids";
- __supportedPrimitives[6] = false;
- __primitivesNames[7] = "second order edges";
- __supportedPrimitives[7] = false;
- __primitivesNames[8] = "second order triangles";
- __supportedPrimitives[8] = false;
- __primitivesNames[9] = "second order quadrangles";
- __supportedPrimitives[9] = false;
- __primitivesNames[10] = "second order tetrahedra";
- __supportedPrimitives[10]= false;
- __primitivesNames[11] = "second order hexahedra";
- __supportedPrimitives[11]= false;
- __primitivesNames[12] = "second order prisms";
- __supportedPrimitives[12]= false;
- __primitivesNames[13] = "second order pyramids";
- __supportedPrimitives[13]= false;
- __primitivesNames[14] = "point";
- __supportedPrimitives[14]= true;
-
- pout() << "Reading file '" << m_filename << "'\n";
-
- // Getting vertices list
- GmshReader::Keyword kw = this->__nextKeyword();
- switch(kw.second) {
- // case NOD: {
- // this->__readGmsh1();
- // break;
- // }
- case MESHFORMAT: {
- double fileVersion = this->_getReal();
- if (fileVersion != 2.2) {
- throw ErrorHandler(__FILE__,__LINE__,
- "Cannot read Gmsh format '"+stringify(fileVersion)+"'",
- ErrorHandler::normal);
- }
- int fileType = this->_getInteger();
- __binary = (fileType == 1);
- if ((fileType < 0) or (fileType > 1)) {
- throw ErrorHandler(__FILE__,__LINE__,
- "Cannot read Gmsh file type '"+stringify(fileType)+"'",
- ErrorHandler::normal);
- }
+ // gmsh 2.2 format keywords
+ __keywordList["$MeshFormat"] = MESHFORMAT;
+ __keywordList["$EndMeshFormat"] = ENDMESHFORMAT;
+ __keywordList["$Nodes"] = NODES;
+ __keywordList["$EndNodes"] = ENDNODES;
+ __keywordList["$Elements"] = ELEMENTS;
+ __keywordList["$EndElements"] = ENDELEMENTS;
+ __keywordList["$PhysicalNames"] = PHYSICALNAMES;
+ __keywordList["$EndPhysicalNames"] = ENDPHYSICALNAMES;
+
+ __numberOfPrimitiveNodes.resize(16);
+ __numberOfPrimitiveNodes[ 0] = 2; // edge
+ __numberOfPrimitiveNodes[ 1] = 3; // triangle
+ __numberOfPrimitiveNodes[ 2] = 4; // quadrangle
+ __numberOfPrimitiveNodes[ 3] = 4; // Tetrahedron
+ __numberOfPrimitiveNodes[ 4] = 8; // Hexaredron
+ __numberOfPrimitiveNodes[ 5] = 6; // Prism
+ __numberOfPrimitiveNodes[ 6] = 5; // Pyramid
+ __numberOfPrimitiveNodes[ 7] = 3; // second order edge
+ __numberOfPrimitiveNodes[ 8] = 6; // second order triangle
+ __numberOfPrimitiveNodes[ 9] = 9; // second order quadrangle
+ __numberOfPrimitiveNodes[10] = 10; // second order tetrahedron
+ __numberOfPrimitiveNodes[11] = 27; // second order hexahedron
+ __numberOfPrimitiveNodes[12] = 18; // second order prism
+ __numberOfPrimitiveNodes[13] = 14; // second order pyramid
+ __numberOfPrimitiveNodes[14] = 1; // point
+
+ __primitivesNames[0] = "edges";
+ __supportedPrimitives[0] = true;
+ __primitivesNames[1] = "triangles";
+ __supportedPrimitives[1] = true;
+ __primitivesNames[2] = "quadrangles";
+ __supportedPrimitives[2] = true;
+ __primitivesNames[3] = "tetrahedra";
+ __supportedPrimitives[3] = true;
+ __primitivesNames[4] = "hexahedra";
+ __supportedPrimitives[4] = true;
+ __primitivesNames[5] = "prisms";
+ __supportedPrimitives[5] = false;
+ __primitivesNames[6] = "pyramids";
+ __supportedPrimitives[6] = false;
+ __primitivesNames[7] = "second order edges";
+ __supportedPrimitives[7] = false;
+ __primitivesNames[8] = "second order triangles";
+ __supportedPrimitives[8] = false;
+ __primitivesNames[9] = "second order quadrangles";
+ __supportedPrimitives[9] = false;
+ __primitivesNames[10] = "second order tetrahedra";
+ __supportedPrimitives[10]= false;
+ __primitivesNames[11] = "second order hexahedra";
+ __supportedPrimitives[11]= false;
+ __primitivesNames[12] = "second order prisms";
+ __supportedPrimitives[12]= false;
+ __primitivesNames[13] = "second order pyramids";
+ __supportedPrimitives[13]= false;
+ __primitivesNames[14] = "point";
+ __supportedPrimitives[14]= true;
+
+ pout() << "Reading file '" << m_filename << "'\n";
+
+ // Getting vertices list
+ GmshReader::Keyword kw = this->__nextKeyword();
+ switch(kw.second) {
+ // case NOD: {
+ // this->__readGmsh1();
+ // break;
+ // }
+ case MESHFORMAT: {
+ double fileVersion = this->_getReal();
+ if (fileVersion != 2.2) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Cannot read Gmsh format '"+stringify(fileVersion)+"'",
+ ErrorHandler::normal);
+ }
+ int fileType = this->_getInteger();
+ __binary = (fileType == 1);
+ if ((fileType < 0) or (fileType > 1)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Cannot read Gmsh file type '"+stringify(fileType)+"'",
+ ErrorHandler::normal);
+ }
- int dataSize = this->_getInteger();
- if (dataSize != sizeof(double)) {
- throw ErrorHandler(__FILE__,__LINE__,
- "Data size not supported '"+stringify(dataSize)+"'",
- ErrorHandler::normal);
- }
+ int dataSize = this->_getInteger();
+ if (dataSize != sizeof(double)) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "Data size not supported '"+stringify(dataSize)+"'",
+ ErrorHandler::normal);
+ }
- if (__binary) {
- // int one=0;
-
- // fseek(__ifh,1L,SEEK_CUR);
- // fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);
-
- // if (one == 1) {
- // __convertEndian = false;
- // } else {
- // invertEndianess(one);
-
- // if (one == 1) {
- // __convertEndian = true;
- // } else {
- // throw ErrorHandler(__FILE__,__LINE__,
- // "Cannot determine endianess",
- // ErrorHandler::normal);
- // }
- // }
-
- // pout() << "- Binary format: ";
- // #ifdef WORDS_BIGENDIAN
- // if (not __convertEndian) {
- // pout() << "big endian\n";
- // } else {
- // pout() << "little endian\n";
- // }
- // #else // WORDS_BIGENDIAN
- // if (not __convertEndian) {
- // pout() << "little endian\n";
- // } else {
- // pout() << "big endian\n";
- // }
- // #endif // WORDS_BIGENDIAN
- }
+ if (__binary) {
+ // int one=0;
+
+ // fseek(__ifh,1L,SEEK_CUR);
+ // fread(reinterpret_cast<char*>(&one),sizeof(int),1,__ifh);
+
+ // if (one == 1) {
+ // __convertEndian = false;
+ // } else {
+ // invertEndianess(one);
+
+ // if (one == 1) {
+ // __convertEndian = true;
+ // } else {
+ // throw ErrorHandler(__FILE__,__LINE__,
+ // "Cannot determine endianess",
+ // ErrorHandler::normal);
+ // }
+ // }
+
+ // pout() << "- Binary format: ";
+ // #ifdef WORDS_BIGENDIAN
+ // if (not __convertEndian) {
+ // pout() << "big endian\n";
+ // } else {
+ // pout() << "little endian\n";
+ // }
+ // #else // WORDS_BIGENDIAN
+ // if (not __convertEndian) {
+ // pout() << "little endian\n";
+ // } else {
+ // pout() << "big endian\n";
+ // }
+ // #endif // WORDS_BIGENDIAN
+ }
- kw = this->__nextKeyword();
- if (kw.second != ENDMESHFORMAT) {
+ kw = this->__nextKeyword();
+ if (kw.second != ENDMESHFORMAT) {
+ throw ErrorHandler(__FILE__,__LINE__,
+ "reading file '"+m_filename
+ +"': expecting $EndMeshFormat, '"+kw.first+"' was found",
+ ErrorHandler::normal);
+ }
+
+ this->__readGmshFormat2_2();
+
+ break;
+ }
+ default: {
throw ErrorHandler(__FILE__,__LINE__,
- "reading file '"+m_filename
- +"': expecting $EndMeshFormat, '"+kw.first+"' was found",
+ "cannot determine format version of '"+m_filename+"'",
ErrorHandler::normal);
}
+ }
- this->__readGmshFormat2_2();
+ this->__proceedData();
+ // this->__createMesh();
+ }
+ catch(const ErrorHandler& e) {
+ e.writeErrorMessage();
+ std::exit(0);
+ }
+ }
- break;
- }
- default: {
- throw ErrorHandler(__FILE__,__LINE__,
- "cannot determine format version of '"+m_filename+"'",
- ErrorHandler::normal);
- }
+ if (commSize() > 1) {
+ pout() << "Sequential mesh read! Need to be dispatched\n";
+ CSRGraph mesh_graph;
+ if (commRank() == 0) {
+ mesh_graph = m_mesh->cellToCellGraph();
}
- this->__proceedData();
- // this->__createMesh();
- }
- catch(const ErrorHandler& e) {
- e.writeErrorMessage();
+ Partitioner P;
+ Array<int> new_cell_owner = P.partition(mesh_graph);
+
+ Messenger::destroy();
std::exit(0);
}
}
diff --git a/src/utils/Partitioner.cpp b/src/utils/Partitioner.cpp
index b3bfd7e705432de9aceb8f5dc8d6a9e34c16ce51..b4c24c5d908019649dcccb52b28ffd81df074091 100644
--- a/src/utils/Partitioner.cpp
+++ b/src/utils/Partitioner.cpp
@@ -2,86 +2,81 @@
#include <Messenger.hpp>
#include <pastis_config.hpp>
+#include <PastisOStream.hpp>
+
#ifdef PASTIS_HAS_MPI
+#define IDXTYPEWIDTH 64
+#define REALTYPEWIDTH 64
#include <parmetis.h>
-#include <vector>
-Partitioner::
-Partitioner()
-{
- std::cout << commRank() << " / " << commSize() << '\n' << std::flush;
- MPI_Barrier(MPI_COMM_WORLD);
- if (commSize() != 3) {
- Messenger::destroy();
- std::cerr << "Commsize=" << commSize() << ". Need 3 PE to test";
- std::exit(1);
- }
+#include <vector>
+
- idx_t wgtflag = 0;
- idx_t numflag = 0;
- idx_t ncon = 1;
- idx_t npart= 3;
- real_t wgts = 1./npart;
- std::vector<real_t> tpwgts{wgts, wgts, wgts};
- std::vector<real_t> ubvec{1.05};
- std::vector<idx_t> options{0,0,0};
- idx_t edgecut = 0;
- std::vector<idx_t> part(5);
+Array<int> Partitioner::partition(const CSRGraph& graph)
+{
+ pout() << "Partitioning graph into "
+ << rang::style::bold << commSize() << rang::style::reset
+ << " parts\n";
MPI_Comm mpi_comm_;
MPI_Comm_dup(MPI_COMM_WORLD, &mpi_comm_);
MPI_Comm* mpi_comm = & mpi_comm_;
- switch (commRank()) {
- case 0: {
- std::vector<idx_t> xadj {0,2,5,8,11,13};
- std::vector<idx_t> adjncy {1,5,0,2,6,1,3,7,2,4,8,3,9};
- std::vector<idx_t> vtxdist{0,5,10,15};
-
- ParMETIS_V3_PartKway(&(vtxdist[0]), &(xadj[0]), &(adjncy[0]),
- NULL, NULL, &wgtflag, &numflag,
- &ncon, &npart, &(tpwgts[0]), &(ubvec[0]),
- &(options[0]), &edgecut, &(part[0]), mpi_comm);
- break;
- }
- case 1: {
- std::vector<idx_t> xadj {0,3,7,11,15,18};
- std::vector<idx_t> adjncy {0,6,10,1,5,7,11,2,6,8,12,3,7,9,13,4,8,14};
- std::vector<idx_t> vtxdist{0,5,10,15};
-
- ParMETIS_V3_PartKway(&(vtxdist[0]), &(xadj[0]), &(adjncy[0]),
- NULL, NULL, &wgtflag, &numflag,
- &ncon, &npart, &(tpwgts[0]), &(ubvec[0]),
- &(options[0]), &edgecut, &(part[0]), mpi_comm);
- break;
- }
- case 2: {
- std::vector<idx_t> xadj {0,2,5,8,11,13};
- std::vector<idx_t> adjncy {5,11,6,10,12,7,11,13,8,12,14,9,13};
- std::vector<idx_t> vtxdist{0,5,10,15};
-
- ParMETIS_V3_PartKway(&(vtxdist[0]), &(xadj[0]), &(adjncy[0]),
- NULL, NULL, &wgtflag, &numflag,
- &ncon, &npart, &(tpwgts[0]), &(ubvec[0]),
- &(options[0]), &edgecut, &(part[0]), mpi_comm);
- break;
- }
- default: {
- std::cerr << "unexpected rank " << commRank() << "!\n";
- std::exit(0);
- }
+
+ int wgtflag = 0;
+ int numflag = 0;
+ int ncon = 1;
+ int npart= commSize();
+ std::vector<float> tpwgts;
+ for (int i_part=0; i_part<npart; ++i_part) {
+ tpwgts.push_back(1./npart);
}
+ std::vector<float> ubvec{1.05};
+ std::vector<int> options{1,1,0};
+ int edgecut = 0;
+ Array<int> part(0);
- std::cerr << commRank() << " FINISHED\n";
- MPI_Barrier(MPI_COMM_WORLD);
+ MPI_Group world_group;
+ MPI_Comm_group(MPI_COMM_WORLD, &world_group);
- Messenger::destroy();
- std::exit(0);
+ MPI_Group mesh_group;
+ std::vector<int> group_ranks{0};
+ MPI_Group_incl(world_group, group_ranks.size(), &(group_ranks[0]), &mesh_group);
+
+ MPI_Comm parmetis_comm;
+ MPI_Comm_create_group(MPI_COMM_WORLD, mesh_group, 1, &parmetis_comm);
+
+ int local_number_of_cells = graph.entries().size()-1;
+
+ if (commRank() ==0) {
+ part = Array<int>(local_number_of_cells);
+ std::vector<int> vtxdist{0,local_number_of_cells};
+
+ static_assert(std::is_same<int, int>());
+
+ const Array<int>& entries = graph.entries();
+ const Array<int>& neighbors = graph.neighbors();
+
+ int result
+ = ParMETIS_V3_PartKway(&(vtxdist[0]), &(entries[0]), &(neighbors[0]),
+ NULL, NULL, &wgtflag, &numflag,
+ &ncon, &npart, &(tpwgts[0]), &(ubvec[0]),
+ &(options[0]), &edgecut, &(part[0]), &parmetis_comm);
+ if (result == METIS_ERROR) {
+ perr() << "Metis Error\n";
+ std::exit(1);
+ }
+ }
+
+ return part;
}
#else // PASTIS_HAS_MPI
-Partitioner::Partitioner() {}
+Array<int> Partitioner::partition(const CSRGraph& graph)
+{
+ return Array<int>(0);
+}
#endif // PASTIS_HAS_MPI
diff --git a/src/utils/Partitioner.hpp b/src/utils/Partitioner.hpp
index 88e3e267250d5a7e8e3b0c37994dc06cca79de1a..427d67f30d4af57f2fe096ffdb8b91736906ee0c 100644
--- a/src/utils/Partitioner.hpp
+++ b/src/utils/Partitioner.hpp
@@ -1,12 +1,16 @@
#ifndef PARTITIONER_HPP
#define PARTITIONER_HPP
+#include <CSRGraph.hpp>
+
class Partitioner
{
public:
- Partitioner();
+ Partitioner() = default;
Partitioner(const Partitioner&) = default;
~Partitioner() = default;
+
+ Array<int> partition(const CSRGraph& graph);
};