From a89d5d1a7f16a483f5e11b9ed0c1d191339d515a Mon Sep 17 00:00:00 2001 From: kyrillh Date: Wed, 17 Apr 2024 10:41:31 +0200 Subject: [PATCH] Structured mesh with overlapping partitioning Small changes --- feddlib/core/FE/Domain_decl.hpp | 2 + feddlib/core/FE/Domain_def.hpp | 5 + feddlib/core/Mesh/MeshPartitioner_decl.hpp | 7 + feddlib/core/Mesh/MeshPartitioner_def.hpp | 1162 ++++++++++------- feddlib/problems/tests/nonLinSchwarz/main.cpp | 66 +- .../tests/nonLinSchwarz/parametersProblem.xml | 2 + 6 files changed, 757 insertions(+), 487 deletions(-) diff --git a/feddlib/core/FE/Domain_decl.hpp b/feddlib/core/FE/Domain_decl.hpp index e901bc6d..3419251e 100644 --- a/feddlib/core/FE/Domain_decl.hpp +++ b/feddlib/core/FE/Domain_decl.hpp @@ -417,6 +417,8 @@ class Domain { */ LO getNumPoints(std::string type="Unique") const;/*local*/ + int getNumProcsCoarseSolve() const;/*local*/ + /*! \brief Checks geometriy @param[in] MeshType diff --git a/feddlib/core/FE/Domain_def.hpp b/feddlib/core/FE/Domain_def.hpp index 060e6126..29f7a2e0 100644 --- a/feddlib/core/FE/Domain_def.hpp +++ b/feddlib/core/FE/Domain_def.hpp @@ -544,6 +544,11 @@ LO Domain::getNumPoints(std::string type) const{ return mesh_->getNumPoints(type); } +template +int Domain::getNumProcsCoarseSolve() const{ + return numProcsCoarseSolve_; +} + template int Domain::checkGeomentry(std::string meshType, int dim) const{ diff --git a/feddlib/core/Mesh/MeshPartitioner_decl.hpp b/feddlib/core/Mesh/MeshPartitioner_decl.hpp index d883bf87..ad186a68 100644 --- a/feddlib/core/Mesh/MeshPartitioner_decl.hpp +++ b/feddlib/core/Mesh/MeshPartitioner_decl.hpp @@ -151,6 +151,13 @@ class MeshPartitioner { */ void readMesh(const int volumeID = 10); + /** + * \brief Build dual graph of the mesh provided where the mesh distributed + * Essentialy a wrapper of the metis function METIS_MeshToDual() + * Note: CSR (compressed sparse row), CRS (compressed row storage), Yale format are all the same + */ + void buildOverlappingDualGraphFromDistributed(const int meshNumber, const int overlap); + /** * \brief Build dual graph of the mesh provided * Essentialy a wrapper of the metis function METIS_MeshToDual() diff --git a/feddlib/core/Mesh/MeshPartitioner_def.hpp b/feddlib/core/Mesh/MeshPartitioner_def.hpp index 37016951..de38e9d1 100644 --- a/feddlib/core/Mesh/MeshPartitioner_def.hpp +++ b/feddlib/core/Mesh/MeshPartitioner_def.hpp @@ -4,6 +4,7 @@ #include "MeshPartitioner_decl.hpp" #include "feddlib/core/FEDDCore.hpp" #include "feddlib/core/LinearAlgebra/Map_decl.hpp" +#include "feddlib/core/LinearAlgebra/MultiVector_decl.hpp" #include "feddlib/core/Utils/FEDDUtils.hpp" #include #include @@ -25,6 +26,8 @@ #include #include #include +#include +#include #include #include #include @@ -36,7 +39,7 @@ /*! Definition of MeshPartitioner - + @brief MeshPartitioner @author Christian Hochmuth @version 1.0 @@ -45,137 +48,130 @@ using namespace std; namespace FEDD { -template -MeshPartitioner::MeshPartitioner() -{ - -} +template MeshPartitioner::MeshPartitioner() {} template -MeshPartitioner::MeshPartitioner( DomainPtrArray_Type domains, ParameterListPtr_Type pL, std::string feType, int dimension ) -{ +MeshPartitioner::MeshPartitioner(DomainPtrArray_Type domains, ParameterListPtr_Type pL, + std::string feType, int dimension) { domains_ = domains; pList_ = pL; feType_ = feType; comm_ = domains_[0]->getComm(); - rankRanges_.resize( domains_.size() ); + rankRanges_.resize(domains_.size()); dim_ = dimension; } - -template -MeshPartitioner::~MeshPartitioner(){ -} +template MeshPartitioner::~MeshPartitioner() {} +template void MeshPartitioner::readAndPartition(int volumeID) { + if (volumeID != 10) { + if (this->comm_->getRank() == 0) { + cout << " #### WARNING: The volumeID was set manually and is no longer 10. Please make sure your volumeID " + "corresponds to the volumeID in your mesh file. #### " + << endl; + } + } + // Read + string delimiter = pList_->get("Delimiter", " "); + for (int i = 0; i < domains_.size(); i++) { + std::string meshName = pList_->get("Mesh " + std::to_string(i + 1) + " Name", "noName"); + TEUCHOS_TEST_FOR_EXCEPTION(meshName == "noName", std::runtime_error, "No mesh name given."); + domains_[i]->initializeUnstructuredMesh(domains_[i]->getDimension(), "P1", + volumeID); // we only allow to read P1 meshes. + domains_[i]->readMeshSize(meshName, delimiter); + } - -template -void MeshPartitioner::readAndPartition( int volumeID) -{ - if(volumeID != 10 ){ - if(this->comm_->getRank()==0){ - cout << " #### WARNING: The volumeID was set manually and is no longer 10. Please make sure your volumeID corresponds to the volumeID in your mesh file. #### " << endl; - } - } - //Read - string delimiter = pList_->get( "Delimiter", " " ); - for (int i=0; iget( "Mesh " + std::to_string(i+1) + " Name", "noName" ); - TEUCHOS_TEST_FOR_EXCEPTION( meshName == "noName", std::runtime_error, "No mesh name given."); - domains_[i]->initializeUnstructuredMesh( domains_[i]->getDimension(), "P1",volumeID ); //we only allow to read P1 meshes. - domains_[i]->readMeshSize( meshName, delimiter ); - } - this->determineRanks(); - for (int i=0; ireadAndPartitionMesh( i ); + for (int i = 0; i < domains_.size(); i++) { + this->readAndPartitionMesh(i); domains_[i]->getMesh()->rankRange_ = rankRanges_[i]; } - } -template -void MeshPartitioner::determineRanks(){ - bool verbose ( comm_->getRank() == 0 ); - vec_int_Type fractions( domains_.size(), 0 ); - bool autoPartition = pList_->get( "Automatic partition", false ); - if ( autoPartition ) { - //determine sum of elements and fractions based of domain contributions +template void MeshPartitioner::determineRanks() { + bool verbose(comm_->getRank() == 0); + vec_int_Type fractions(domains_.size(), 0); + bool autoPartition = pList_->get("Automatic partition", false); + if (autoPartition) { + // determine sum of elements and fractions based of domain contributions GO sumElements = 0; - for (int i=0; igetNumElementsGlobal(); - - for (int i=0; igetNumElementsGlobal()*100) / sumElements; + + for (int i = 0; i < fractions.size(); i++) + fractions[i] = (domains_[i]->getNumElementsGlobal() * 100) / sumElements; int diff = std::accumulate(fractions.begin(), fractions.end(), 0) - 100; auto iterator = fractions.begin(); - while (diff>0) { + while (diff > 0) { (*iterator)--; iterator++; diff--; } iterator = fractions.begin(); - while (diff<0) { + while (diff < 0) { (*iterator)++; iterator++; diff++; } - this->determineRanksFromFractions( fractions ); - + this->determineRanksFromFractions(fractions); + if (verbose) { std::cout << "\t --- ---------------- ---" << std::endl; std::cout << "\t --- Mesh Partitioner ---" << std::endl; std::cout << "\t --- ---------------- ---" << std::endl; - std::cout << "\t --- Automatic partition for "<< comm_->getSize() <<" ranks" << std::endl; - for (int i=0; i( rankRanges_[i] )<< " to "<< get<1>( rankRanges_[i] ) << std::endl; + std::cout << "\t --- Automatic partition for " << comm_->getSize() << " ranks" << std::endl; + for (int i = 0; i < domains_.size(); i++) { + std::cout << "\t --- Fraction mesh " << to_string(i + 1) << " : " << fractions[i] + << " of 100; rank range: " << get<0>(rankRanges_[i]) << " to " << get<1>(rankRanges_[i]) + << std::endl; } } - - } - else if( autoPartition == false && pList_->get("Mesh 1 fraction ranks",-1) >= 0 ){ - for (int i=0; fractions.size(); i++) - fractions[i] = pList_->get("Mesh " + std::to_string(i+1) + " fraction ranks", -1); - - TEUCHOS_TEST_FOR_EXCEPTION( std::accumulate(fractions.begin(), fractions.end(), 0) != 100, std::runtime_error, "Fractions do not sum up to 100!"); - this->determineRanksFromFractions( fractions ); - + + } else if (autoPartition == false && pList_->get("Mesh 1 fraction ranks", -1) >= 0) { + for (int i = 0; fractions.size(); i++) + fractions[i] = pList_->get("Mesh " + std::to_string(i + 1) + " fraction ranks", -1); + + TEUCHOS_TEST_FOR_EXCEPTION(std::accumulate(fractions.begin(), fractions.end(), 0) != 100, std::runtime_error, + "Fractions do not sum up to 100!"); + this->determineRanksFromFractions(fractions); + if (verbose) { std::cout << "\t --- ---------------- ---" << std::endl; std::cout << "\t --- Mesh Partitioner ---" << std::endl; std::cout << "\t --- ---------------- ---" << std::endl; - std::cout << "\t --- Fraction partition for "<< comm_->getSize() <<" ranks" << std::endl; - for (int i=0; i( rankRanges_[i] )<< " to "<< get<1>( rankRanges_[i] ) << std::endl; + std::cout << "\t --- Fraction partition for " << comm_->getSize() << " ranks" << std::endl; + for (int i = 0; i < domains_.size(); i++) { + std::cout << "\t --- Fraction mesh " << to_string(i + 1) << " : " << fractions[i] + << " of 100; rank range: " << get<0>(rankRanges_[i]) << " to " << get<1>(rankRanges_[i]) + << std::endl; } } - } - else if( autoPartition == false && pList_->get("Mesh 1 fraction ranks",-1) < 0 && pList_->get("Mesh 1 number ranks",-1) > 0 ){ + } else if (autoPartition == false && pList_->get("Mesh 1 fraction ranks", -1) < 0 && + pList_->get("Mesh 1 number ranks", -1) > 0) { int size = comm_->getSize(); vec_int_Type numberRanks(domains_.size()); - for (int i=0; iget("Mesh " + std::to_string(i+1) + " number ranks",0); - - TEUCHOS_TEST_FOR_EXCEPTION( std::accumulate(numberRanks.begin(), numberRanks.end(), 0) > size, std::runtime_error, "Too many ranks requested for mesh partition!"); - this->determineRanksFromNumberRanks( numberRanks ); + for (int i = 0; i < numberRanks.size(); i++) + numberRanks[i] = pList_->get("Mesh " + std::to_string(i + 1) + " number ranks", 0); + + TEUCHOS_TEST_FOR_EXCEPTION(std::accumulate(numberRanks.begin(), numberRanks.end(), 0) > size, + std::runtime_error, "Too many ranks requested for mesh partition!"); + this->determineRanksFromNumberRanks(numberRanks); if (verbose) { std::cout << "\t --- ---------------- ---" << std::endl; std::cout << "\t --- Mesh Partitioner ---" << std::endl; std::cout << "\t --- ---------------- ---" << std::endl; - std::cout << "\t --- Rank number partition for "<< comm_->getSize() <<" ranks" << std::endl; - for (int i=0; i( rankRanges_[i] )<< " to "<< get<1>( rankRanges_[i] ) << std::endl; + std::cout << "\t --- Rank number partition for " << comm_->getSize() << " ranks" << std::endl; + for (int i = 0; i < domains_.size(); i++) { + std::cout << "\t --- Rank range mesh " << to_string(i + 1) << " :" << get<0>(rankRanges_[i]) << " to " + << get<1>(rankRanges_[i]) << std::endl; } } - - } - else{ - for (int i=0; igetSize()-1 ); + + } else { + for (int i = 0; i < domains_.size(); i++) + rankRanges_[i] = std::make_tuple(0, comm_->getSize() - 1); if (verbose) { std::cout << "\t --- ---------------- ---" << std::endl; std::cout << "\t --- Mesh Partitioner ---" << std::endl; @@ -186,207 +182,204 @@ void MeshPartitioner::determineRanks(){ } template -void MeshPartitioner::determineRanksFromFractions( vec_int_Type& fractions ){ - - int lowerRank = 0; int size = comm_->getSize(); +void MeshPartitioner::determineRanksFromFractions(vec_int_Type &fractions) { + + int lowerRank = 0; + int size = comm_->getSize(); int upperRank = 0; - for (int i=0; i1) + rankRanges_[i] = std::make_tuple(lowerRank, upperRank); + if (size > 1) lowerRank = upperRank + 1; else lowerRank = upperRank; } int startLoc = 0; - while ( upperRank > size-1 ) { - for (int i=startLoc; i0) - std::get<0>( rankRanges_[i] )--; - std::get<1>( rankRanges_[i] )--; + while (upperRank > size - 1) { + for (int i = startLoc; i < rankRanges_.size(); i++) { + if (i > 0) + std::get<0>(rankRanges_[i])--; + std::get<1>(rankRanges_[i])--; } startLoc++; upperRank--; } startLoc = 0; - while ( upperRank < size-1 ) { - for (int i=startLoc; i0) - std::get<0>( rankRanges_[i] )++; - std::get<1>( rankRanges_[i] )++; + while (upperRank < size - 1) { + for (int i = startLoc; i < rankRanges_.size(); i++) { + if (i > 0) + std::get<0>(rankRanges_[i])++; + std::get<1>(rankRanges_[i])++; } startLoc++; upperRank++; } } - + template -void MeshPartitioner::determineRanksFromNumberRanks( vec_int_Type& numberRanks ){ - - int lowerRank = 0; int size = comm_->getSize(); +void MeshPartitioner::determineRanksFromNumberRanks(vec_int_Type &numberRanks) { + + int lowerRank = 0; + int size = comm_->getSize(); int upperRank = 0; - for (int i=0; i size-1 ) { - for (int i=startLoc; i0) - std::get<0>( rankRanges_[i] )--; - std::get<1>( rankRanges_[i] )--; + while (upperRank > size - 1) { + for (int i = startLoc; i < rankRanges_.size(); i++) { + if (i > 0) + std::get<0>(rankRanges_[i])--; + std::get<1>(rankRanges_[i])--; } startLoc++; upperRank--; } startLoc = 0; - while ( upperRank < size-1 ) { - for (int i=startLoc; i0) - std::get<0>( rankRanges_[i] )++; - std::get<1>( rankRanges_[i] )++; + while (upperRank < size - 1) { + for (int i = startLoc; i < rankRanges_.size(); i++) { + if (i > 0) + std::get<0>(rankRanges_[i])++; + std::get<1>(rankRanges_[i])++; } startLoc++; upperRank++; } - } -/// Reading and partioning of the mesh. Input File is .mesh. Reading is serial and at some point the mesh entities are distributed along the processors. +/// Reading and partioning of the mesh. Input File is .mesh. Reading is serial and at some point the mesh entities are +/// distributed along the processors. template -void MeshPartitioner::readAndPartitionMesh( int meshNumber ){ - +void MeshPartitioner::readAndPartitionMesh(int meshNumber) { + typedef Teuchos::OrdinalTraits OTGO; #ifdef UNDERLYING_LIB_TPETRA string underlyingLib = "Tpetra"; #endif - - MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast( domains_[meshNumber]->getMesh() ); - - // Reading nodes + + MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast(domains_[meshNumber]->getMesh()); + + // Reading nodes meshUnstr->readMeshEntity("node"); - // We delete the point at this point. We only need the flags to determine surface elements. We will load them again later. + // We delete the point at this point. We only need the flags to determine surface elements. We will load them again + // later. meshUnstr->pointsRep_.reset(); // Reading elements meshUnstr->readMeshEntity("element"); // Reading surfaces meshUnstr->readMeshEntity("surface"); - // Reading line segments + // Reading line segments meshUnstr->readMeshEntity("line"); - - - bool verbose ( comm_->getRank() == 0 ); + bool verbose(comm_->getRank() == 0); bool buildEdges = pList_->get("Build Edge List", true); bool buildSurfaces = pList_->get("Build Surface List", true); - // Adding surface as subelement to elements + // Adding surface as subelement to elements if (buildSurfaces) - this->setSurfacesToElements( meshNumber ); + this->setSurfacesToElements(meshNumber); else meshUnstr->deleteSurfaceElements(); - - // Serially distributed elements + + // Serially distributed elements ElementsPtr_Type elementsMesh = meshUnstr->getElementsC(); - + // Setup Metis idx_t options[METIS_NOPTIONS]; METIS_SetDefaultOptions(options); options[METIS_OPTION_NUMBERING] = 0; options[METIS_OPTION_SEED] = 666; - options[METIS_OPTION_CONTIG] = pList_->get("Contiguous",false); //0: Does not force contiguous partitions; 1: Forces contiguous partitions. - options[METIS_OPTION_MINCONN] = 0; // 1: Explicitly minimize the maximum connectivity. - options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT;// or METIS_OBJTYPE_VOL + options[METIS_OPTION_CONTIG] = + pList_->get("Contiguous", false); // 0: Does not force contiguous partitions; 1: Forces contiguous partitions. + options[METIS_OPTION_MINCONN] = 0; // 1: Explicitly minimize the maximum connectivity. + options[METIS_OPTION_OBJTYPE] = METIS_OBJTYPE_CUT; // or METIS_OBJTYPE_VOL // options[METIS_OPTION_RTYPE] = METIS_RTYPE_GREEDY; options[METIS_OPTION_NITER] = 50; // default is 10 options[METIS_OPTION_CCORDER] = 1; - idx_t ne = meshUnstr->getNumElementsGlobal(); // Global number of elements - idx_t nn = meshUnstr->getNumGlobalNodes(); // Global number of nodes + idx_t ne = meshUnstr->getNumElementsGlobal(); // Global number of elements + idx_t nn = meshUnstr->getNumGlobalNodes(); // Global number of nodes idx_t ned = meshUnstr->getEdgeElements()->numberElements(); // Global number of edges - int dim = meshUnstr->getDimension(); std::string FEType = domains_[meshNumber]->getFEType(); - // Setup for paritioning with metis + // Setup for paritioning with metis vec_idx_Type eptr_vec(0); // Vector for local elements ptr (local is still global at this point) vec_idx_Type eind_vec(0); // Vector for local elements ids - + makeContinuousElements(elementsMesh, eind_vec, eptr_vec); idx_t *eptr = &eptr_vec.at(0); idx_t *eind = &eind_vec.at(0); - + idx_t ncommon; int orderSurface; - if (dim==2) { - if (FEType=="P1") { + if (dim == 2) { + if (FEType == "P1") { ncommon = 2; - } - else if(FEType=="P2"){ + } else if (FEType == "P2") { ncommon = 3; } - } - else if (dim==3) { - if (FEType=="P1") { + } else if (dim == 3) { + if (FEType == "P1") { ncommon = 3; - } - else if(FEType=="P2"){ + } else if (FEType == "P2") { ncommon = 6; } - } - else + } else TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong Dimension."); - + idx_t objval = 0; - vec_idx_Type epart(ne,-1); - vec_idx_Type npart(nn,-1); + vec_idx_Type epart(ne, -1); + vec_idx_Type npart(nn, -1); // Partitioning elements with metis if (verbose) cout << "-- Start partitioning with Metis ... " << flush; - + { - FEDD_TIMER_START(partitionTimer," : MeshPartitioner : Partition Elements"); - idx_t nparts = std::get<1>( rankRanges_[meshNumber] ) - std::get<0>( rankRanges_[meshNumber] ) + 1; - if ( nparts > 1 ) { + FEDD_TIMER_START(partitionTimer, " : MeshPartitioner : Partition Elements"); + idx_t nparts = std::get<1>(rankRanges_[meshNumber]) - std::get<0>(rankRanges_[meshNumber]) + 1; + if (nparts > 1) { int rank = this->comm_->getRank(); // upperRange - lowerRange +1 - idx_t returnCode = METIS_PartMeshDual(&ne, &nn, eptr, eind, NULL, NULL, &ncommon, &nparts, NULL, options, &objval, &epart[0], &npart[0]); - if ( verbose ) + idx_t returnCode = METIS_PartMeshDual(&ne, &nn, eptr, eind, NULL, NULL, &ncommon, &nparts, NULL, options, + &objval, &epart[0], &npart[0]); + if (verbose) cout << "\n--\t Metis return code: " << returnCode; - } - else{ - for (int i=0; igetRank() - std::get<0>( rankRanges_[meshNumber] ) ){ + // Getting global IDs of element's nodes + for (int i = 0; i < ne; i++) { + if (epart[i] == comm_->getRank() - std::get<0>(rankRanges_[meshNumber])) { locepart.push_back(i); - for (int j=eptr[i]; j::readAndPartitionMesh( int meshNumber ){ make_unique(pointsRepIndices); if (verbose) cout << "done!" << endl; - + // Building repeated node map - Teuchos::ArrayView pointsRepGlobMapping = Teuchos::arrayViewFromVector( pointsRepIndices ); - meshUnstr->mapRepeated_.reset( new Map( underlyingLib, OTGO::invalid(), pointsRepGlobMapping, 0, this->comm_) ); + Teuchos::ArrayView pointsRepGlobMapping = Teuchos::arrayViewFromVector(pointsRepIndices); + meshUnstr->mapRepeated_.reset( + new Map(underlyingLib, OTGO::invalid(), pointsRepGlobMapping, 0, this->comm_)); MapConstPtr_Type mapRepeated = meshUnstr->mapRepeated_; // We keep the global elements if we want to build edges later. Otherwise they will be deleted - ElementsPtr_Type elementsGlobal = Teuchos::rcp( new Elements_Type( *elementsMesh ) ); + ElementsPtr_Type elementsGlobal = Teuchos::rcp(new Elements_Type(*elementsMesh)); - // Resetting elements to add the corresponding local IDs instead of global ones - meshUnstr->elementsC_.reset(new Elements ( FEType, dim ) ); + // Resetting elements to add the corresponding local IDs instead of global ones + meshUnstr->elementsC_.reset(new Elements(FEType, dim)); { - Teuchos::ArrayView elementsGlobalMapping = Teuchos::arrayViewFromVector( locepart ); + Teuchos::ArrayView elementsGlobalMapping = Teuchos::arrayViewFromVector(locepart); // elementsGlobalMapping -> elements per Processor - meshUnstr->elementMap_.reset(new Map( underlyingLib, (GO) -1, elementsGlobalMapping, 0, this->comm_) ); - + meshUnstr->elementMap_.reset(new Map(underlyingLib, (GO)-1, elementsGlobalMapping, 0, this->comm_)); + { int localSurfaceCounter = 0; - for (int i=0; i tmpElement; - for (int j=eptr[locepart.at(i)]; jgetLocalElement( (long long) eind[j] ); + for (int j = eptr[locepart.at(i)]; j < eptr[locepart.at(i) + 1]; j++) { + // local indices + int index = mapRepeated->getLocalElement((long long)eind[j]); tmpElement.push_back(index); } - //std::sort(tmpElement.begin(), tmpElement.end()); - FiniteElement fe( tmpElement, elementsGlobal->getElement( locepart.at(i) ).getFlag() ); + // std::sort(tmpElement.begin(), tmpElement.end()); + FiniteElement fe(tmpElement, elementsGlobal->getElement(locepart.at(i)).getFlag()); // convert global IDs of (old) globally owned subelements to local IDs if (buildSurfaces) { - FiniteElement feGlobalIDs = elementsGlobal->getElement( locepart.at(i) ); - if (feGlobalIDs.subElementsInitialized()){ + FiniteElement feGlobalIDs = elementsGlobal->getElement(locepart.at(i)); + if (feGlobalIDs.subElementsInitialized()) { ElementsPtr_Type subEl = feGlobalIDs.getSubElements(); - subEl->globalToLocalIDs( mapRepeated ); - fe.setSubElements( subEl ); + subEl->globalToLocalIDs(mapRepeated); + fe.setSubElements(subEl); } } - meshUnstr->elementsC_->addElement( fe ); + meshUnstr->elementsC_->addElement(fe); } } } - // Next we distribute the coordinates and flags correctly - - + // Next we distribute the coordinates and flags correctly + meshUnstr->readMeshEntity("node"); // We reread the nodes, as they were deleted earlier - + if (verbose) cout << "-- Build Repeated Points Volume ... " << flush; - + // building the unique map - meshUnstr->mapUnique_ = meshUnstr->mapRepeated_->buildUniqueMap( rankRanges_[meshNumber] ); + meshUnstr->mapUnique_ = meshUnstr->mapRepeated_->buildUniqueMap(rankRanges_[meshNumber]); // free(epart); if (verbose) @@ -456,166 +449,166 @@ void MeshPartitioner::readAndPartitionMesh( int meshNumber ){ { vec2D_dbl_Type points = *meshUnstr->getPointsRepeated(); vec_int_Type flags = *meshUnstr->getBCFlagRepeated(); - meshUnstr->pointsRep_.reset( new std::vector >( meshUnstr->mapRepeated_->getNodeNumElements(), std::vector(dim,-1.) ) ); - meshUnstr->bcFlagRep_.reset( new std::vector ( meshUnstr->mapRepeated_->getNodeNumElements(), 0 ) ); - + meshUnstr->pointsRep_.reset(new std::vector>(meshUnstr->mapRepeated_->getNodeNumElements(), + std::vector(dim, -1.))); + meshUnstr->bcFlagRep_.reset(new std::vector(meshUnstr->mapRepeated_->getNodeNumElements(), 0)); + int pointIDcont; - for (int i=0; ipointsRep_->at(i).at(j) = points[pointIDcont][j]; - meshUnstr->bcFlagRep_->at(i) = flags[pointIDcont]; + meshUnstr->bcFlagRep_->at(i) = flags[pointIDcont]; } } - // Setting unique points and flags - meshUnstr->pointsUni_.reset(new std::vector >( meshUnstr->mapUnique_->getNodeNumElements(), std::vector(dim,-1.) ) ); - meshUnstr->bcFlagUni_.reset(new std::vector ( meshUnstr->mapUnique_->getNodeNumElements(), 0) ); + // Setting unique points and flags + meshUnstr->pointsUni_.reset(new std::vector>(meshUnstr->mapUnique_->getNodeNumElements(), + std::vector(dim, -1.))); + meshUnstr->bcFlagUni_.reset(new std::vector(meshUnstr->mapUnique_->getNodeNumElements(), 0)); GO indexGlobal; MapConstPtr_Type map = meshUnstr->getMapRepeated(); vec2D_dbl_ptr_Type pointsRep = meshUnstr->pointsRep_; - for (int i=0; imapUnique_->getNodeNumElements() ; i++) { + for (int i = 0; i < meshUnstr->mapUnique_->getNodeNumElements(); i++) { indexGlobal = meshUnstr->mapUnique_->getGlobalElement(i); - for (int j=0; jpointsUni_->at(i).at(j) = pointsRep->at( map->getLocalElement( indexGlobal) ).at(j); + for (int j = 0; j < dim; j++) { + meshUnstr->pointsUni_->at(i).at(j) = pointsRep->at(map->getLocalElement(indexGlobal)).at(j); } - meshUnstr->bcFlagUni_->at(i) = meshUnstr->bcFlagRep_->at( map->getLocalElement( indexGlobal) ); + meshUnstr->bcFlagUni_->at(i) = meshUnstr->bcFlagRep_->at(map->getLocalElement(indexGlobal)); } - // Finally we build the edges. As part of the edge build involves nodes and elements, - // they should be build last to avoid any local and global IDs mix up - if (!buildEdges) + // Finally we build the edges. As part of the edge build involves nodes and elements, + // they should be build last to avoid any local and global IDs mix up + if (!buildEdges) elementsGlobal.reset(); - locepart.erase(locepart.begin(),locepart.end()); + locepart.erase(locepart.begin(), locepart.end()); if (verbose) cout << "done!" << endl; - - if (buildSurfaces){ - this->setEdgesToSurfaces( meshNumber ); // Adding edges as subelements in the 3D case. All dim-1-Subelements were already set - } - else + + if (buildSurfaces) { + this->setEdgesToSurfaces( + meshNumber); // Adding edges as subelements in the 3D case. All dim-1-Subelements were already set + } else meshUnstr->deleteSurfaceElements(); - + if (buildEdges) { if (verbose) cout << "-- Build edge element list ... \n" << flush; - - buildEdgeListParallel( meshUnstr, elementsGlobal ); - + + buildEdgeListParallel(meshUnstr, elementsGlobal); + if (verbose) - cout << "\n done!"<< endl; - - MapConstPtr_Type elementMap = meshUnstr->getElementMap(); + cout << "\n done!" << endl; + + MapConstPtr_Type elementMap = meshUnstr->getElementMap(); - FEDD_TIMER_START(partitionEdgesTimer," : MeshPartitioner : Partition Edges"); - meshUnstr->getEdgeElements()->partitionEdges( elementMap, mapRepeated ); + FEDD_TIMER_START(partitionEdgesTimer, " : MeshPartitioner : Partition Edges"); + meshUnstr->getEdgeElements()->partitionEdges(elementMap, mapRepeated); FEDD_TIMER_STOP(partitionEdgesTimer); // edge global indices on different processors - for( int i=0; igetEdgeElements()->numberElements() ; i++){ - locedpart.push_back(meshUnstr->getEdgeElements()->getGlobalID((LO) i)); + for (int i = 0; i < meshUnstr->getEdgeElements()->numberElements(); i++) { + locedpart.push_back(meshUnstr->getEdgeElements()->getGlobalID((LO)i)); } // Setup for the EdgeMap - Teuchos::ArrayView edgesGlobalMapping = Teuchos::arrayViewFromVector( locedpart ); - meshUnstr->edgeMap_.reset(new Map( underlyingLib, (GO) -1, edgesGlobalMapping, 0, this->comm_) ); + Teuchos::ArrayView edgesGlobalMapping = Teuchos::arrayViewFromVector(locedpart); + meshUnstr->edgeMap_.reset(new Map(underlyingLib, (GO)-1, edgesGlobalMapping, 0, this->comm_)); } - if (verbose) cout << "done!" << endl; - + if (verbose) cout << "-- Partition interface ... " << flush; meshUnstr->partitionInterface(); - + if (verbose) cout << "done!" << endl; - } - -/// Function that (addionally) adds edges as subelements in the 3D case. This is relevant when the .mesh file also contains edge information -/// in 3D. Otherwise edges are not set as subelements in 3D. +/// Function that (addionally) adds edges as subelements in the 3D case. This is relevant when the .mesh file also +/// contains edge information in 3D. Otherwise edges are not set as subelements in 3D. template -void MeshPartitioner::setEdgesToSurfaces(int meshNumber){ - bool verbose ( comm_->getRank() == 0 ); - MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast( domains_[meshNumber]->getMesh() ); +void MeshPartitioner::setEdgesToSurfaces(int meshNumber) { + bool verbose(comm_->getRank() == 0); + MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast(domains_[meshNumber]->getMesh()); ElementsPtr_Type elementsMesh = meshUnstr->getElementsC(); MapConstPtr_Type mapRepeated = meshUnstr->mapRepeated_; if (verbose) cout << "-- Set edges of surfaces of elements ... " << flush; - - FEDD_TIMER_START(surfacesTimer," : MeshPartitioner : Set Surfaces of Edge Elements"); + + FEDD_TIMER_START(surfacesTimer, " : MeshPartitioner : Set Surfaces of Edge Elements"); vec2D_int_Type localEdgeIDPermutation; - setLocalSurfaceEdgeIndices( localEdgeIDPermutation, meshUnstr->getEdgeElementOrder() ); - + setLocalSurfaceEdgeIndices(localEdgeIDPermutation, meshUnstr->getEdgeElementOrder()); + int volumeID = meshUnstr->volumeID_; - + ElementsPtr_Type elements = meshUnstr->getElementsC(); ElementsPtr_Type edgeElements = meshUnstr->getSurfaceEdgeElements(); - + /* First, we convert the surface edge elements to a 2D array so we can use std::find.*/ // Can we use/implement find for Elements_Type? - vec2D_int_Type edgeElements_vec( edgeElements->numberElements() ); - vec_int_Type edgeElementsFlag_vec( edgeElements->numberElements() ); - for (int i=0; inumberElements()); + vec_int_Type edgeElementsFlag_vec(edgeElements->numberElements()); + for (int i = 0; i < edgeElements_vec.size(); i++) { vec_int_Type edge = edgeElements->getElement(i).getVectorNodeListNonConst(); - std::sort( edge.begin(), edge.end() ); - edgeElements_vec.at(i) = edge; + std::sort(edge.begin(), edge.end()); + edgeElements_vec.at(i) = edge; edgeElementsFlag_vec.at(i) = edgeElements->getElement(i).getFlag(); } - - vec_int_ptr_Type flags = meshUnstr->bcFlagRep_; + + vec_int_ptr_Type flags = meshUnstr->bcFlagRep_; int elementEdgeSurfaceCounter; - for (int i=0; inumberElements(); i++) { + for (int i = 0; i < elements->numberElements(); i++) { elementEdgeSurfaceCounter = 0; bool mark = false; - for (int j=0; jgetElement( i ).size(); j++) { - if ( flags->at( elements->getElement( i ).getNode( j ) ) < volumeID ) + for (int j = 0; j < elements->getElement(i).size(); j++) { + if (flags->at(elements->getElement(i).getNode(j)) < volumeID) elementEdgeSurfaceCounter++; } - if (elementEdgeSurfaceCounter >= meshUnstr->getEdgeElementOrder()){ - //We want to find all surfaces of element i and set the surfaces to the element - findAndSetSurfaceEdges( edgeElements_vec, edgeElementsFlag_vec, elements->getElement(i), localEdgeIDPermutation, mapRepeated ); + if (elementEdgeSurfaceCounter >= meshUnstr->getEdgeElementOrder()) { + // We want to find all surfaces of element i and set the surfaces to the element + findAndSetSurfaceEdges(edgeElements_vec, edgeElementsFlag_vec, elements->getElement(i), + localEdgeIDPermutation, mapRepeated); } } if (verbose) cout << "done!" << endl; } -/// Adding surfaces as subelements to the corresponding elements. This adresses surfaces in 3D or edges in 2D. -/// If in the 3D case edges are also part of the .mesh file, they will be added as subelements later. +/// Adding surfaces as subelements to the corresponding elements. This adresses surfaces in 3D or edges in 2D. +/// If in the 3D case edges are also part of the .mesh file, they will be added as subelements later. template -void MeshPartitioner::setSurfacesToElements(int meshNumber){ - - bool verbose ( comm_->getRank() == 0 ); - MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast( domains_[meshNumber]->getMesh() ); - ElementsPtr_Type elementsMesh = meshUnstr->getElementsC(); // Previously read Elements +void MeshPartitioner::setSurfacesToElements(int meshNumber) { + + bool verbose(comm_->getRank() == 0); + MeshUnstrPtr_Type meshUnstr = Teuchos::rcp_dynamic_cast(domains_[meshNumber]->getMesh()); + ElementsPtr_Type elementsMesh = meshUnstr->getElementsC(); // Previously read Elements if (verbose) - cout << "-- Set surfaces of elements ... " << flush; - - FEDD_TIMER_START(surfacesTimer," : MeshPartitioner : Set Surfaces of Elements"); - + cout << "-- Set surfaces of elements ... " << flush; + + FEDD_TIMER_START(surfacesTimer, " : MeshPartitioner : Set Surfaces of Elements"); + vec2D_int_Type localSurfaceIDPermutation; // get permutations - setLocalSurfaceIndices( localSurfaceIDPermutation, meshUnstr->getSurfaceElementOrder() ); - + setLocalSurfaceIndices(localSurfaceIDPermutation, meshUnstr->getSurfaceElementOrder()); + int volumeID = meshUnstr->volumeID_; ElementsPtr_Type surfaceElements = meshUnstr->getSurfaceElements(); - + /* First, we convert the surface Elements to a 2D array so we can use std::find. and partition it at the same time. We use a simple linear partition. This is done to reduce times spend for std::find. The result is then communicated. We therefore use the unpartitoned elements*/ // Can we use/implement find for Elements_Type? - + int size = this->comm_->getSize(); - - LO numSurfaceEl = surfaceElements->numberElements();// / size; + + LO numSurfaceEl = surfaceElements->numberElements(); // / size; /*LO rest = surfaceElements->numberElements() % size; - + vec_GO_Type offsetVec(size); for (int i=0; i::setSurfacesToElements(int meshNumber){ offsetVec[i]+=rest; }*/ - vec2D_int_Type surfElements_vec( numSurfaceEl ); - vec2D_int_Type surfElements_vec_sorted( numSurfaceEl ); + vec2D_int_Type surfElements_vec(numSurfaceEl); + vec2D_int_Type surfElements_vec_sorted(numSurfaceEl); - vec_int_Type surfElementsFlag_vec( numSurfaceEl ); + vec_int_Type surfElementsFlag_vec(numSurfaceEl); vec_GO_Type offsetVec(size); int offset = offsetVec[this->comm_->getRank()]; - for (int i=0; igetElement(i ).getVectorNodeListNonConst(); // surfaceElements->getElement(i + offset).getVectorNodeListNonConst(); - surfElements_vec.at(i) = surface; - std::sort( surface.begin(), surface.end() ); // We need to maintain a consistent numbering in the surface elements, so we use a sorted and unsorted vector + for (int i = 0; i < surfElements_vec.size(); i++) { + vec_int_Type surface = + surfaceElements->getElement(i) + .getVectorNodeListNonConst(); // surfaceElements->getElement(i + offset).getVectorNodeListNonConst(); + surfElements_vec.at(i) = surface; + std::sort(surface.begin(), surface.end()); // We need to maintain a consistent numbering in the surface + // elements, so we use a sorted and unsorted vector surfElements_vec_sorted.at(i) = surface; - surfElementsFlag_vec.at(i) = surfaceElements->getElement(i).getFlag(); // surfaceElements->getElement(i + offset).getFlag(); - + surfElementsFlag_vec.at(i) = + surfaceElements->getElement(i).getFlag(); // surfaceElements->getElement(i + offset).getFlag(); } - + // Delete the surface elements. They will be added to the elements in the following loop. surfaceElements.reset(); - vec_int_ptr_Type flags = meshUnstr->bcFlagRep_; - + vec_int_ptr_Type flags = meshUnstr->bcFlagRep_; + int elementSurfaceCounter; int surfaceElOrder = meshUnstr->getSurfaceElementOrder(); - for (int i=0; inumberElements(); i++) { + for (int i = 0; i < elementsMesh->numberElements(); i++) { elementSurfaceCounter = 0; - for (int j=0; jgetElement( i ).size(); j++) { - if ( flags->at( elementsMesh->getElement( i ).getNode( j ) ) < volumeID ) - elementSurfaceCounter++; + for (int j = 0; j < elementsMesh->getElement(i).size(); j++) { + if (flags->at(elementsMesh->getElement(i).getNode(j)) < volumeID) + elementSurfaceCounter++; } - - if (elementSurfaceCounter >= surfaceElOrder){ - FEDD_TIMER_START(findSurfacesTimer," : MeshPartitioner : Find and Set Surfaces"); - //We want to find all surfaces of element i and set the surfaces to the element - this->findAndSetSurfacesPartitioned( surfElements_vec_sorted, surfElements_vec, surfElementsFlag_vec, elementsMesh->getElement(i), localSurfaceIDPermutation, offsetVec, i ); + + if (elementSurfaceCounter >= surfaceElOrder) { + FEDD_TIMER_START(findSurfacesTimer, " : MeshPartitioner : Find and Set Surfaces"); + // We want to find all surfaces of element i and set the surfaces to the element + this->findAndSetSurfacesPartitioned(surfElements_vec_sorted, surfElements_vec, surfElementsFlag_vec, + elementsMesh->getElement(i), localSurfaceIDPermutation, offsetVec, i); } } - + if (verbose) cout << "done!" << endl; } -/// Function that sets edges (2D) or surfaces(3D) to the corresponding element. It determines for the specific element 'element' if it has a surface with a non -/// volumeflag that can be set as subelement +/// Function that sets edges (2D) or surfaces(3D) to the corresponding element. It determines for the specific element +/// 'element' if it has a surface with a non volumeflag that can be set as subelement template -void MeshPartitioner::findAndSetSurfacesPartitioned( vec2D_int_Type& surfElements_vec, vec2D_int_Type& surfElements_vec_unsorted, vec_int_Type& surfElementsFlag_vec, FiniteElement& element, vec2D_int_Type& permutation , vec_GO_Type& linearSurfacePartitionOffset, int globalElID ){ - - // In general we look through the different permutations the input element 'element' can have and if they correspond to a surface. - // The mesh's surface elements 'surfElements_vec' are then used to determine the corresponding surface - // If found, the nodes are then used to build the subelement and the corresponding surfElementFlag is set. - // The Ids are global at this point, as the values are not distributed yet. +void MeshPartitioner::findAndSetSurfacesPartitioned( + vec2D_int_Type &surfElements_vec, vec2D_int_Type &surfElements_vec_unsorted, vec_int_Type &surfElementsFlag_vec, + FiniteElement &element, vec2D_int_Type &permutation, vec_GO_Type &linearSurfacePartitionOffset, int globalElID) { + + // In general we look through the different permutations the input element 'element' can have and if they correspond + // to a surface. The mesh's surface elements 'surfElements_vec' are then used to determine the corresponding surface + // If found, the nodes are then used to build the subelement and the corresponding surfElementFlag is set. + // The Ids are global at this point, as the values are not distributed yet. int loc, id1Glob, id2Glob, id3Glob; int size = this->comm_->getSize(); vec_int_Type locAll(size); if (dim_ == 2) { - for (int j=0; j id2Glob){ + if (id1Glob > id2Glob) { tmpSurface[0] = id2Glob; tmpSurface[1] = id1Glob; - } - else{ + } else { tmpSurface[0] = id1Glob; tmpSurface[1] = id2Glob; } - - loc = searchInSurfaces( surfElements_vec, tmpSurface ); - - Teuchos::gatherAll( *this->comm_, 1, &loc, locAll.size(), &locAll[0] ); + + loc = searchInSurfaces(surfElements_vec, tmpSurface); + + Teuchos::gatherAll(*this->comm_, 1, &loc, locAll.size(), &locAll[0]); int surfaceRank = -1; int counter = 0; - while (surfaceRank<0 && counter -1) surfaceRank = counter; counter++; } int surfFlag = -1; - if (loc>-1) + if (loc > -1) surfFlag = surfElementsFlag_vec[loc]; - - if (surfaceRank>-1) { - Teuchos::broadcast(*this->comm_,surfaceRank,1,&loc); - Teuchos::broadcast(*this->comm_,surfaceRank,1,&surfFlag); - - FiniteElement feSurface( tmpSurface, surfFlag); - if ( !element.subElementsInitialized() ) - element.initializeSubElements( "P1", 1 ); // only P1 for now - - element.addSubElement( feSurface ); + + if (surfaceRank > -1) { + Teuchos::broadcast(*this->comm_, surfaceRank, 1, &loc); + Teuchos::broadcast(*this->comm_, surfaceRank, 1, &surfFlag); + + FiniteElement feSurface(tmpSurface, surfFlag); + if (!element.subElementsInitialized()) + element.initializeSubElements("P1", 1); // only P1 for now + + element.addSubElement(feSurface); } } - } - else if (dim_ == 3){ - for (int j=0; j( *this->comm_, 1, &loc, locAll.size(), &locAll[0] ); - + int surfaceRank = -1; int counter = 0; while (surfaceRank<0 && counter::findAndSetSurfacesPartitioned( vec2D_int_Type int surfFlag = -1; if (loc>-1) surfFlag = surfElementsFlag_vec[loc]; - + if (surfaceRank>-1) {*/ - if(loc > -1 ){ - //Teuchos::broadcast(*this->comm_,surfaceRank,1,&loc); - //Teuchos::broadcast(*this->comm_,surfaceRank,1,&surfFlag); + if (loc > -1) { + // Teuchos::broadcast(*this->comm_,surfaceRank,1,&loc); + // Teuchos::broadcast(*this->comm_,surfaceRank,1,&surfFlag); int surfFlag = surfElementsFlag_vec[loc]; - //cout << " Surfaces set to elements on Proc " << this->comm_->getRank() << " " << surfElements_vec_unsorted[loc][0] << " " << surfElements_vec_unsorted[loc][1] << " " << surfElements_vec_unsorted[loc][2] << endl; - FiniteElement feSurface( surfElements_vec_unsorted[loc], surfFlag); - if ( !element.subElementsInitialized() ) - element.initializeSubElements( "P1", 2 ); // only P1 for now - - element.addSubElement( feSurface ); + // cout << " Surfaces set to elements on Proc " << this->comm_->getRank() << " " << + // surfElements_vec_unsorted[loc][0] << " " << surfElements_vec_unsorted[loc][1] << " " << + // surfElements_vec_unsorted[loc][2] << endl; + FiniteElement feSurface(surfElements_vec_unsorted[loc], surfFlag); + if (!element.subElementsInitialized()) + element.initializeSubElements("P1", 2); // only P1 for now + + element.addSubElement(feSurface); } } } - } - template -void MeshPartitioner::buildEdgeListParallel( MeshUnstrPtr_Type mesh, ElementsPtr_Type elementsGlobal ){ - FEDD_TIMER_START(edgeListTimer," : MeshReader : Build Edge List"); +void MeshPartitioner::buildEdgeListParallel(MeshUnstrPtr_Type mesh, ElementsPtr_Type elementsGlobal) { + FEDD_TIMER_START(edgeListTimer, " : MeshReader : Build Edge List"); ElementsPtr_Type elements = mesh->getElementsC(); - - TEUCHOS_TEST_FOR_EXCEPTION( elements->getFiniteElementType() != "P1", std::runtime_error ,"Unknown discretization for method buildEdgeList(...)."); + + TEUCHOS_TEST_FOR_EXCEPTION(elements->getFiniteElementType() != "P1", std::runtime_error, + "Unknown discretization for method buildEdgeList(...)."); CommConstPtr_Type comm = mesh->getComm(); - bool verbose(comm->getRank()==0); - - MapConstPtr_Type repeatedMap = mesh->getMapRepeated(); + bool verbose(comm->getRank() == 0); + + MapConstPtr_Type repeatedMap = mesh->getMapRepeated(); // Building local edges with repeated node list vec2D_int_Type localEdgeIndices(0); - setLocalEdgeIndices( localEdgeIndices ); - EdgeElementsPtr_Type edges = Teuchos::rcp( new EdgeElements_Type() ); - for (int i=0; inumberElements(); i++) { - for (int j=0; jgetElement( i ).getNode( localEdgeIndices[j][0] ); - int id2 = elementsGlobal->getElement( i ).getNode( localEdgeIndices[j][1] ); + setLocalEdgeIndices(localEdgeIndices); + EdgeElementsPtr_Type edges = Teuchos::rcp(new EdgeElements_Type()); + for (int i = 0; i < elementsGlobal->numberElements(); i++) { + for (int j = 0; j < localEdgeIndices.size(); j++) { + + int id1 = elementsGlobal->getElement(i).getNode(localEdgeIndices[j][0]); + int id2 = elementsGlobal->getElement(i).getNode(localEdgeIndices[j][1]); vec_int_Type edgeVec(2); - if (id1addEdge( edge, i ); + + FiniteElement edge(edgeVec); + edges->addEdge(edge, i); } } // we do not need elementsGlobal anymore elementsGlobal.reset(); - + vec2D_GO_Type combinedEdgeElements; - FEDD_TIMER_START(edgeListUniqueTimer," : MeshReader : Make Edge List Unique"); - edges->sortUniqueAndSetGlobalIDs( combinedEdgeElements ); + FEDD_TIMER_START(edgeListUniqueTimer, " : MeshReader : Make Edge List Unique"); + edges->sortUniqueAndSetGlobalIDs(combinedEdgeElements); FEDD_TIMER_STOP(edgeListUniqueTimer); - //Next we need to communicate all edge information. This will not scale at all! - - edges->setElementsEdges( combinedEdgeElements ); - - mesh->setEdgeElements( edges ); - + // Next we need to communicate all edge information. This will not scale at all! + + edges->setElementsEdges(combinedEdgeElements); + + mesh->setEdgeElements(edges); }; template -void MeshPartitioner::setLocalEdgeIndices(vec2D_int_Type &localEdgeIndices ){ - if ( dim_ == 2 ) { - localEdgeIndices.resize(3,vec_int_Type(2,-1)); +void MeshPartitioner::setLocalEdgeIndices(vec2D_int_Type &localEdgeIndices) { + if (dim_ == 2) { + localEdgeIndices.resize(3, vec_int_Type(2, -1)); localEdgeIndices.at(0).at(0) = 0; localEdgeIndices.at(0).at(1) = 1; localEdgeIndices.at(1).at(0) = 0; localEdgeIndices.at(1).at(1) = 2; localEdgeIndices.at(2).at(0) = 1; localEdgeIndices.at(2).at(1) = 2; - } - else if( dim_ == 3) { - localEdgeIndices.resize(6,vec_int_Type(2,-1)); + } else if (dim_ == 3) { + localEdgeIndices.resize(6, vec_int_Type(2, -1)); localEdgeIndices.at(0).at(0) = 0; localEdgeIndices.at(0).at(1) = 1; localEdgeIndices.at(1).at(0) = 0; @@ -837,35 +832,34 @@ void MeshPartitioner::setLocalEdgeIndices(vec2D_int_Type &localEdge localEdgeIndices.at(4).at(1) = 3; localEdgeIndices.at(5).at(0) = 2; localEdgeIndices.at(5).at(1) = 3; - } } - + template -void MeshPartitioner::makeContinuousElements(ElementsPtr_Type elements, vec_idx_Type& eind_vec, vec_idx_Type& eptr_vec ){ - +void MeshPartitioner::makeContinuousElements(ElementsPtr_Type elements, vec_idx_Type &eind_vec, + vec_idx_Type &eptr_vec) { + int nodesPerElement = elements->nodesPerElement(); - int elcounter=0; - for (int i=0; inumberElements(); i++) { - for (int j=0; jgetElement( i ).getNode( j ) ); + int elcounter = 0; + for (int i = 0; i < elements->numberElements(); i++) { + for (int j = 0; j < nodesPerElement; j++) { + eind_vec.push_back(elements->getElement(i).getNode(j)); } eptr_vec.push_back(elcounter); elcounter += nodesPerElement; } eptr_vec.push_back(elcounter); - } - template -void MeshPartitioner::setLocalSurfaceEdgeIndices( vec2D_int_Type &localSurfaceEdgeIndices, int edgesElementOrder ){ - - if ( dim_ == 3 ) { - - if (edgesElementOrder == 2) { //P1 - localSurfaceEdgeIndices.resize(6,vec_int_Type(2,-1)); +void MeshPartitioner::setLocalSurfaceEdgeIndices(vec2D_int_Type &localSurfaceEdgeIndices, + int edgesElementOrder) { + + if (dim_ == 3) { + + if (edgesElementOrder == 2) { // P1 + localSurfaceEdgeIndices.resize(6, vec_int_Type(2, -1)); localSurfaceEdgeIndices.at(0).at(0) = 0; localSurfaceEdgeIndices.at(0).at(1) = 1; localSurfaceEdgeIndices.at(1).at(0) = 0; @@ -883,89 +877,91 @@ void MeshPartitioner::setLocalSurfaceEdgeIndices( vec2D_int_Type &l } template -void MeshPartitioner::findAndSetSurfaceEdges( vec2D_int_Type& edgeElements_vec, vec_int_Type& edgeElementsFlag_vec, FiniteElement& element, vec2D_int_Type& permutation, MapConstPtr_Type mapRepeated){ - +void MeshPartitioner::findAndSetSurfaceEdges(vec2D_int_Type &edgeElements_vec, + vec_int_Type &edgeElementsFlag_vec, FiniteElement &element, + vec2D_int_Type &permutation, + MapConstPtr_Type mapRepeated) { + int loc, id1Glob, id2Glob; - if (dim_ == 3){ - for (int j=0; jgetGlobalElement( element.getNode( permutation.at(j).at(0) ) ); - id2Glob = mapRepeated->getGlobalElement( element.getNode( permutation.at(j).at(1) ) ); + if (dim_ == 3) { + for (int j = 0; j < permutation.size(); j++) { + id1Glob = mapRepeated->getGlobalElement(element.getNode(permutation.at(j).at(0))); + id2Glob = mapRepeated->getGlobalElement(element.getNode(permutation.at(j).at(1))); vec_int_Type tmpEdge(0); if (id2Glob > id1Glob) - tmpEdge = {id1Glob , id2Glob}; + tmpEdge = {id1Glob, id2Glob}; else - tmpEdge = {id2Glob , id1Glob}; - - loc = searchInSurfaces( edgeElements_vec, tmpEdge ); - - if (loc>-1) { - - int id1 = element.getNode( permutation.at(j).at(0) ); - int id2 = element.getNode( permutation.at(j).at(1) ); + tmpEdge = {id2Glob, id1Glob}; + + loc = searchInSurfaces(edgeElements_vec, tmpEdge); + + if (loc > -1) { + + int id1 = element.getNode(permutation.at(j).at(0)); + int id2 = element.getNode(permutation.at(j).at(1)); vec_int_Type tmpEdgeLocal(0); - if (id2>id1) - tmpEdgeLocal = { id1 , id2 }; + if (id2 > id1) + tmpEdgeLocal = {id1, id2}; else - tmpEdgeLocal = { id2 , id1 }; - - // If no partition was performed, all information is still global at this point. We still use the function below and partition the mesh and surfaces later. - FiniteElement feEdge( tmpEdgeLocal, edgeElementsFlag_vec[loc] ); - // In some cases an edge is the only part of the surface of an Element. In that case there does not exist a triangle subelement. - // We then have to initialize the edge as subelement. - // In very coarse meshes it is even possible that an interior element has multiple surface edges or nodes connected to the surface, which is why we might even set interior edges as subelements - if ( !element.subElementsInitialized() ){ - element.initializeSubElements( "P1", 1 ); // only P1 for now - element.addSubElement( feEdge ); - } - else{ + tmpEdgeLocal = {id2, id1}; + + // If no partition was performed, all information is still global at this point. We still use the + // function below and partition the mesh and surfaces later. + FiniteElement feEdge(tmpEdgeLocal, edgeElementsFlag_vec[loc]); + // In some cases an edge is the only part of the surface of an Element. In that case there does not + // exist a triangle subelement. We then have to initialize the edge as subelement. In very coarse meshes + // it is even possible that an interior element has multiple surface edges or nodes connected to the + // surface, which is why we might even set interior edges as subelements + if (!element.subElementsInitialized()) { + element.initializeSubElements("P1", 1); // only P1 for now + element.addSubElement(feEdge); + } else { ElementsPtr_Type surfaces = element.getSubElements(); - if(surfaces->getDimension() == 2) // We set the edge to the corresponding surface element(s) - surfaces->setToCorrectElement( feEdge ); // Case we have surface subelements - else{ // simply adding the edge as subelement - element.addSubElement( feEdge ); // Case we dont have surface subelements - // Comment: Theoretically edges as subelement are only truely relevant when we apply a neuman bc on an edge which is currently not the case. - // This is just in case! + if (surfaces->getDimension() == 2) // We set the edge to the corresponding surface element(s) + surfaces->setToCorrectElement(feEdge); // Case we have surface subelements + else { // simply adding the edge as subelement + element.addSubElement(feEdge); // Case we dont have surface subelements + // Comment: Theoretically edges as subelement are only truely relevant when we apply a neuman bc + // on an edge which is currently not the case. This is just in case! } - } + } } } } - } template -int MeshPartitioner::searchInSurfaces( vec2D_int_Type& surfaces, vec_int_Type searchSurface){ - +int MeshPartitioner::searchInSurfaces(vec2D_int_Type &surfaces, vec_int_Type searchSurface) { + int loc = -1; - - vec2D_int_Type::iterator it = find(surfaces.begin(),surfaces.end(), searchSurface); - - if (it!=surfaces.end()) - loc = distance(surfaces.begin(),it); - + + vec2D_int_Type::iterator it = find(surfaces.begin(), surfaces.end(), searchSurface); + + if (it != surfaces.end()) + loc = distance(surfaces.begin(), it); + return loc; } - + template -void MeshPartitioner::setLocalSurfaceIndices(vec2D_int_Type &localSurfaceIndices, int surfaceElementOrder ){ - - if ( dim_ == 2 ) { - - if (surfaceElementOrder == 2) { //P1 - localSurfaceIndices.resize(3,vec_int_Type(3,-1)); +void MeshPartitioner::setLocalSurfaceIndices(vec2D_int_Type &localSurfaceIndices, + int surfaceElementOrder) { + + if (dim_ == 2) { + + if (surfaceElementOrder == 2) { // P1 + localSurfaceIndices.resize(3, vec_int_Type(3, -1)); localSurfaceIndices.at(0).at(0) = 0; localSurfaceIndices.at(0).at(1) = 1; localSurfaceIndices.at(1).at(0) = 0; localSurfaceIndices.at(1).at(1) = 2; localSurfaceIndices.at(2).at(0) = 1; localSurfaceIndices.at(2).at(1) = 2; - } - else + } else TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No permutation for this surface yet."); - } - else if ( dim_ == 3 ){ + } else if (dim_ == 3) { if (surfaceElementOrder == 3) { - localSurfaceIndices.resize(4,vec_int_Type(3,-1)); + localSurfaceIndices.resize(4, vec_int_Type(3, -1)); localSurfaceIndices.at(0).at(0) = 0; localSurfaceIndices.at(0).at(1) = 1; localSurfaceIndices.at(0).at(2) = 2; @@ -978,8 +974,7 @@ void MeshPartitioner::setLocalSurfaceIndices(vec2D_int_Type &localS localSurfaceIndices.at(3).at(0) = 0; localSurfaceIndices.at(3).at(1) = 2; localSurfaceIndices.at(3).at(2) = 3; - } - else + } else TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "No permutation for this surface yet."); } } @@ -1025,6 +1020,223 @@ template void MeshPartitioner +void MeshPartitioner::buildOverlappingDualGraphFromDistributed(const int meshNumber, const int overlap) { + + typedef Teuchos::OrdinalTraits OTGO; + const auto myRank = this->comm_->getRank(); + +#ifdef UNDERLYING_LIB_TPETRA + const Xpetra::UnderlyingLib underlyingLibType = Xpetra::UseTpetra; + const string underlyingLib = "Tpetra"; +#endif + const auto mesh = this->domains_[meshNumber]->getMesh(); + TEUCHOS_TEST_FOR_EXCEPTION(mesh.is_null(), std::runtime_error, "No mesh to work with."); + + // Number of elements in the mesh + auto ne = Teuchos::implicit_cast(mesh->getElementMap()->getGlobalNumElements()); + // Number of nodes in the mesh + auto nn = Teuchos::implicit_cast(mesh->getMapUnique()->getGlobalNumElements()); + + // Setup for paritioning with metis + vec_idx_Type eptrVec(0); // Vector for local elements ptr (local is still global at this point) + vec_idx_Type eindVec(0); // Vector for local elements ids + // Serially distributed elements + const auto elementsMesh = mesh->getElementsC(); + // Fill the arrays + makeContinuousElements(elementsMesh, eindVec, eptrVec); + + // Make the arrays locally replicated for METIS + // Build distributed vector for element indices. They live on the elementMap. + auto eptrVecDistributed = Teuchos::rcp(new MultiVector(mesh->getElementMap())); + // Ignore the last entry (size) so that existing element map can be used + for (auto i = 0; i < eptrVec.size() - 1; i++) { + eptrVecDistributed->replaceLocalValue(i, 0, eptrVec.at(i) + myRank * eptrVec.back()); + } + + // Build distributed vector for node indices. They live on the repeated map. + auto eindMap = Teuchos::rcp(new Map( + underlyingLib, eindVec.size() * (comm_->getSize() - this->domains_[0]->getNumProcsCoarseSolve()), + eindVec.size(), 0, comm_)); + auto eindVecDistributed = Teuchos::rcp(new MultiVector(eindMap, 1)); + + for (auto i = 0; i < eindVec.size(); i++) { + eindVecDistributed->replaceLocalValue(i, 0, mesh->getMapRepeated()->getGlobalElement(eindVec.at(i))); + } + + // Build locally replicated maps + auto locReplNodeMap = Teuchos::rcp(new Map(Xpetra::MapFactory::createLocalMap( + underlyingLibType, eindVecDistributed->getXpetraMultiVector()->getGlobalLength(), this->comm_))); + auto locReplElemMap = Teuchos::rcp( + new Map(Xpetra::MapFactory::createLocalMap(underlyingLibType, ne, this->comm_))); + + auto eptrVecLocRepl = Teuchos::rcp(new MultiVector(locReplElemMap, 1)); + eptrVecLocRepl->importFromVector(eptrVecDistributed); + auto eindVecLocRepl = Teuchos::rcp(new MultiVector(locReplNodeMap, 1)); + eindVecLocRepl->importFromVector(eindVecDistributed); + + eindVec.resize(eindVecLocRepl->getLocalLength()); + for (auto i = 0; i < eindVec.size(); i++) { + eindVec.at(i) = eindVecLocRepl->getData(0)[i]; + } + + eptrVec.resize(eptrVecLocRepl->getLocalLength() + 1); + for (auto i = 0; i < eptrVec.size(); i++) { + eptrVec.at(i) = eptrVecLocRepl->getData(0)[i]; + } + // Add the missing size entry + eptrVec.back() = eindVec.size(); + + auto eptr = &eptrVec.at(0); + auto eind = &eindVec.at(0); + + // Number of common nodes required to classify two elements as neighbours: min(ncommon, n1-1, n2-1) + idx_t ncommon; + const auto dim = mesh->getDimension(); + const auto FEType = domains_[meshNumber]->getFEType(); + if (dim == 2) { + if (FEType == "P1") { + ncommon = 2; + } else if (FEType == "P2") { + ncommon = 3; + } + } else if (dim == 3) { + if (FEType == "P1") { + ncommon = 3; + } else if (FEType == "P2") { + ncommon = 6; + } + } else { + TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong Dimension."); + } + idx_t numflag = 0; + + // Arrays for storing the partition + // METIS allocates these using malloc + idx_t *xadj; + idx_t *adjncy; + + if (myRank == 0) { + cout << "--- Building dual graph with METIS ..."; + } + const auto returnCode = METIS_MeshToDual(&ne, &nn, eptr, eind, &ncommon, &numflag, &xadj, &adjncy); + + if (myRank == 0) { + cout << "\n--\t Metis return code: " << returnCode; + cout << "\n--- done" << endl; + } + + /* mesh->elementMap_.reset(new Map(locReplElemMap->getXpetraMap())); */ + + // Copy idx_t arrays to ArrayRCP arrays required by graph constructor + auto xadjArrayRCP = Teuchos::ArrayRCP(ne + 1); + for (auto i = 0; i < xadjArrayRCP.size(); i++) { + xadjArrayRCP[i] = xadj[i]; + } + // Size of adjncy is stored in last entry of xadj + auto adjncyArrayRCP = Teuchos::ArrayRCP(xadj[ne]); + for (auto i = 0; i < adjncyArrayRCP.size(); i++) { + adjncyArrayRCP[i] = adjncy[i]; + } + + // Both row and column maps are unity + // All entries are on all ranks + mesh->dualGraph_ = Xpetra::CrsGraphFactory::Build( + locReplElemMap->getXpetraMap(), locReplElemMap->getXpetraMap(), xadjArrayRCP, adjncyArrayRCP); + mesh->dualGraph_->fillComplete(); + METIS_Free(xadj); + METIS_Free(adjncy); + + // ########### Partition the dual graph and extend overlap as specified ############### + // We need the distributed dual graph to be able to perform assembly on each subdomain + + // Create an export object since the target map is one to one but the source map may not be + auto exporter = + Xpetra::ExportFactory::Build(locReplElemMap->getXpetraMap(), mesh->getElementMap()->getXpetraMap()); + + // New graph into which the previous graph will be imported + auto newDualGraph = Xpetra::CrsGraphFactory::Build( + mesh->getElementMap()->getXpetraMap(), mesh->getElementMap()->getXpetraMap()->getLocalNumElements()); + + newDualGraph->doExport(*mesh->dualGraph_, *exporter, Xpetra::INSERT); + newDualGraph->fillComplete(); + + // Cast to pointer-to-const for FROSch function + auto graphExtended = Teuchos::rcp_implicit_cast>(newDualGraph); + + // Extend overlap by specified number of times + for (auto i = 0; i < overlap; i++) { + ExtendOverlapByOneLayer(graphExtended, graphExtended); + } + + // Store the interior of the subdomain to differentiate the border + auto extendedElementMap = FROSch::SortMapByGlobalIndex(graphExtended->getRowMap()); + + // Build graph with sorted map + newDualGraph = + Xpetra::CrsGraphFactory::Build(extendedElementMap, extendedElementMap->getLocalNumElements()); + auto importer = Xpetra::ImportFactory::Build(graphExtended->getRowMap(), extendedElementMap); + newDualGraph->doImport(*graphExtended, *importer, Xpetra::ADD); + newDualGraph->fillComplete(graphExtended->getDomainMap(), graphExtended->getRangeMap()); + + mesh->dualGraph_ = newDualGraph; + + // Rebuild elements to be locally replicated. Very inefficient approach since they will be rebuilt on the + // overlapping map later but makes use of existing functions + auto flagsDistributed = Teuchos::rcp(new MultiVector(mesh->getElementMap(), 1)); + for (auto i = 0; i < mesh->getElementMap()->getNodeNumElements(); i++) { + flagsDistributed->getDataNonConst(0)[i] = elementsMesh->getElement(i).getFlag(); + } + auto flagsLocRepl = Teuchos::rcp(new MultiVector(locReplElemMap, 1)); + flagsLocRepl->importFromVector(flagsDistributed); + + mesh->elementsC_.reset(new Elements(FEType, dim)); + for (auto i = 0; i < locReplElemMap->getNodeNumElements(); i++) { + // Build local element and save + std::vector tmpElement; + // Iterate over the nodes in an element + for (int j = eptrVec.at(i); j < eptrVec.at(i + 1); j++) { + // Map the node index from global to local and save + tmpElement.push_back(eindVec.at(j)); + } + FiniteElement tempFE(tmpElement, flagsLocRepl->getData(0)[i]); + // NOTE KHo Surfaces are not added here for now. Since they probably will not be needed. + mesh->elementsC_->addElement(tempFE); + } + // Communicate points and store in points repeated + auto pointsRepDistributed = Teuchos::rcp(new MultiVector(mesh->getMapRepeated(), dim)); + for (auto i = 0; i < mesh->pointsRep_->size(); i++) { + for (auto j = 0; j < dim; j++) { + pointsRepDistributed->getDataNonConst(j)[i] = mesh->pointsRep_->at(i).at(j); + } + } + auto locReplPointMap = Teuchos::rcp(new Map(Xpetra::MapFactory::createLocalMap( + underlyingLibType, mesh->getMapUnique()->getGlobalNumElements(), this->comm_))); + + auto pointsRepLocRepl = Teuchos::rcp(new MultiVector(locReplPointMap, dim)); + pointsRepLocRepl->importFromVector(pointsRepDistributed); + mesh->pointsRep_.reset( + new std::vector>(locReplPointMap->getNodeNumElements(), std::vector(dim, -1.))); + + for (auto i = 0; i < mesh->pointsRep_->size(); i++) { + for (auto j = 0; j < dim; j++) { + mesh->pointsRep_->at(i).at(j) = pointsRepLocRepl->getData(j)[i]; + } + } + // Communicate boundary condition flags and store in BCFlagRepeated + auto bcFlagsDistributed = Teuchos::rcp(new MultiVector(mesh->getMapRepeated(), 1)); + for (auto i = 0; i < mesh->bcFlagRep_->size(); i++) { + bcFlagsDistributed->getDataNonConst(0)[i] = mesh->bcFlagRep_->at(i); + } + auto bcFlagsRepLocRepl = Teuchos::rcp(new MultiVector(locReplPointMap, 1)); + bcFlagsRepLocRepl->importFromVector(bcFlagsDistributed); + mesh->bcFlagRep_.reset(new std::vector(locReplPointMap->getNodeNumElements(), 0)); + + for (auto i = 0; i < mesh->bcFlagRep_->size(); i++) { + mesh->bcFlagRep_->at(i) = bcFlagsRepLocRepl->getData(0)[i]; + } +} + template void MeshPartitioner::buildDualGraph(const int meshNumber) { @@ -1279,13 +1491,13 @@ void MeshPartitioner::buildSubdomainFEsAndNodeLists(const int me const Xpetra::UnderlyingLib underlyingLibType = Xpetra::UseTpetra; #endif - const auto myRank = comm_->getRank(); - const auto meshUnstr = Teuchos::rcp_dynamic_cast(domains_.at(meshNumber)->getMesh()); - const auto dim = meshUnstr->getDimension(); - const auto FEType = domains_.at(meshNumber)->getFEType(); - const auto elementMap = meshUnstr->getElementMap(); - const auto elementMapOverlapping = meshUnstr->dualGraph_->getRowMap(); - const int indexBase = 0; + auto myRank = comm_->getRank(); + auto meshUnstr = domains_.at(meshNumber)->getMesh(); + auto dim = meshUnstr->getDimension(); + auto FEType = domains_.at(meshNumber)->getFEType(); + auto elementMap = meshUnstr->getElementMap(); + auto elementMapOverlapping = meshUnstr->dualGraph_->getRowMap(); + int indexBase = 0; vec_int_Type elementsOverlappingGhostsIndices(0); // Build index vectors for the elements which can then be used to build the elements in the current subdomain. diff --git a/feddlib/problems/tests/nonLinSchwarz/main.cpp b/feddlib/problems/tests/nonLinSchwarz/main.cpp index 41f9f61c..6bc8efdc 100644 --- a/feddlib/problems/tests/nonLinSchwarz/main.cpp +++ b/feddlib/problems/tests/nonLinSchwarz/main.cpp @@ -2,6 +2,7 @@ #include "feddlib/core/General/DefaultTypeDefs.hpp" #include "feddlib/core/LinearAlgebra/BlockMultiVector_decl.hpp" #include "feddlib/core/Mesh/MeshPartitioner_decl.hpp" +#include "feddlib/core/Utils/FEDDUtils.hpp" #include "feddlib/problems/Solver/NonLinearSolver.hpp" #include "feddlib/problems/specific/NonLinLaplace_decl.hpp" #include @@ -69,11 +70,27 @@ int main(int argc, char *argv[]) { // ######################## // Set default values for command line parameters // ######################## + // Command Line Parameters + Teuchos::CommandLineProcessor myCLP; + string ulib_str = "Tpetra"; + myCLP.setOption("ulib",&ulib_str,"Underlying lib"); - bool debug = false; - if (argc > 1) { - debug = true; + string xmlProblemFile = "parametersProblem.xml"; + myCLP.setOption("problemfile",&xmlProblemFile,".xml file with Inputparameters."); + string xmlPrecFile = "parametersPrec.xml"; + myCLP.setOption("precfile",&xmlPrecFile,".xml file with Inputparameters."); + string xmlSolverFile = "parametersSolver.xml"; + myCLP.setOption("solverfile",&xmlSolverFile,".xml file with Inputparameters."); + + myCLP.recogniseAllOptions(true); + myCLP.throwExceptions(false); + Teuchos::CommandLineProcessor::EParseCommandLineReturn parseReturn = myCLP.parse(argc,argv); + if(parseReturn == Teuchos::CommandLineProcessor::PARSE_HELP_PRINTED) { + mpiSession.~GlobalMPISession(); + return 0; } + + bool debug = false; if (comm->getRank() == 1 && debug) { waitForGdbAttach(); } @@ -87,10 +104,6 @@ int main(int argc, char *argv[]) { cout << "###############################################################" << endl; } - string xmlProblemFile = "parametersProblem.xml"; - string xmlPrecFile = "parametersPrec.xml"; - string xmlSolverFile = "parametersSolver.xml"; - ParameterListPtr_Type parameterListProblem = Teuchos::getParametersFromXmlFile(xmlProblemFile); ParameterListPtr_Type parameterListPrec = Teuchos::getParametersFromXmlFile(xmlPrecFile); ParameterListPtr_Type parameterListSolver = Teuchos::getParametersFromXmlFile(xmlSolverFile); @@ -100,8 +113,14 @@ int main(int argc, char *argv[]) { parameterListAll->setParameters(*parameterListSolver); int dim = parameterListProblem->sublist("Parameter").get("Dimension", 2); + int minNumberSubdomains = 1; + int n = -1; + int m = parameterListProblem->sublist("Parameter").get("H/h", 5); + string meshType = parameterListProblem->sublist("Parameter").get("Mesh Type", "structured"); string FEType = parameterListProblem->sublist("Parameter").get("Discretization", "P1"); auto overlap = parameterListSolver->get("Overlap", 1); + int numProcsCoarseSolve = 0; // TODO: kho need to implement this + int size = comm->getSize() - numProcsCoarseSolve; // ######################## // Build mesh @@ -116,11 +135,35 @@ int main(int argc, char *argv[]) { ParameterListPtr_Type pListPartitioner = sublist(parameterListAll, "Mesh Partitioner"); MeshPartitioner partitionerP1(domainP1Array, pListPartitioner, "P1", dim); - partitionerP1.readMesh(); - partitionerP1.buildDualGraph(0); - partitionerP1.partitionDualGraphWithOverlap(0, overlap); + if (!meshType.compare("structured")) { + TEUCHOS_TEST_FOR_EXCEPTION(size % minNumberSubdomains != 0, std::logic_error, + "Wrong number of processors for structured mesh."); + if (dim == 2) { + n = (int)(std::pow(size, 1 / 2.) + 100. * Teuchos::ScalarTraits::eps()); // 1/H + std::vector x(2); + x[0] = 0.0; + x[1] = 0.0; + domainP1.reset(new Domain(x, 1., 1., comm)); + domainP1Array[0] = domainP1; + domainP1->buildMesh(1, "Square", dim, FEType, n, m, numProcsCoarseSolve); + } else if (dim == 3) { + n = (int)(std::pow(size, 1 / 3.) + 100. * Teuchos::ScalarTraits::eps()); // 1/H + std::vector x(3); + x[0] = 0.0; + x[1] = 0.0; + x[2] = 0.0; + domainP1.reset(new Domain(x, 1., 1., 1., comm)); + domainP1Array[0] = domainP1; + domainP1->buildMesh(1, "Square", dim, FEType, n, m, numProcsCoarseSolve); + } + partitionerP1 = MeshPartitioner(domainP1Array, pListPartitioner, "P1", dim); + partitionerP1.buildOverlappingDualGraphFromDistributed(0, 2); + } else if (!meshType.compare("unstructured")) { + partitionerP1.readMesh(); + partitionerP1.buildDualGraph(0); + partitionerP1.partitionDualGraphWithOverlap(0, overlap); + } partitionerP1.buildSubdomainFEsAndNodeLists(0); - domain = domainP1; // ######################## @@ -158,7 +201,6 @@ int main(int argc, char *argv[]) { nlSolverAssFE.solve(*nonLinLaplace); comm->barrier(); - // Export Solution bool boolExportSolution = true; if (boolExportSolution) { diff --git a/feddlib/problems/tests/nonLinSchwarz/parametersProblem.xml b/feddlib/problems/tests/nonLinSchwarz/parametersProblem.xml index 3151657e..a1b5c09c 100644 --- a/feddlib/problems/tests/nonLinSchwarz/parametersProblem.xml +++ b/feddlib/problems/tests/nonLinSchwarz/parametersProblem.xml @@ -1,6 +1,8 @@ + +