46#ifndef MUELU_LOCALLEXICOGRAPHICINDEXMANAGER_DEF_HPP_
47#define MUELU_LOCALLEXICOGRAPHICINDEXMANAGER_DEF_HPP_
50#include <Xpetra_MapFactory.hpp>
54 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
57 const int NumDimensions,
const int interpolationOrder,
58 const int MyRank,
const int NumRanks,
59 const Array<GO> GFineNodesPerDir,
const Array<LO> LFineNodesPerDir,
60 const Array<LO> CoarseRate,
const Array<GO> MeshData) :
61 IndexManager(comm, coupled, false, NumDimensions, interpolationOrder, GFineNodesPerDir, LFineNodesPerDir),
62 myRank(MyRank), numRanks(NumRanks) {
70 for(
int dim = 0; dim < 3; ++dim) {
72 if(CoarseRate.size() == 1) {
74 }
else if(CoarseRate.size() == this->numDimensions) {
83 for(
int rank = 0; rank <
numRanks; ++rank) {
85 for(
int entry = 0; entry < 10; ++entry) {
86 meshData[rank][entry] = MeshData[10*rank + entry];
97 for(
int dim = 0; dim < 3; ++dim) {
107 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 this->gNumCoarseNodes10 = this->gCoarseNodesPerDir[0]*this->gCoarseNodesPerDir[1];
111 this->gNumCoarseNodes = this->gNumCoarseNodes10*this->gCoarseNodesPerDir[2];
114 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 Array<LO>& ghostedNodeCoarseLIDs,
118 Array<int>& ghostedNodeCoarsePIDs,
119 Array<GO>& ghostedNodeCoarseGIDs)
const {
122 ghostedNodeCoarseLIDs.resize(this->getNumLocalGhostedNodes());
123 ghostedNodeCoarsePIDs.resize(this->getNumLocalGhostedNodes());
124 ghostedNodeCoarseGIDs.resize(this->numGhostedNodes);
131 Array<LO> ghostedCoarseNodeCoarseIndices(3), ghostedCoarseNodeFineIndices(3);
132 Array<LO> lCoarseNodeCoarseIndices(3);
133 Array<GO> lCoarseNodeCoarseGIDs(this->lNumCoarseNodes);
134 LO currentIndex = -1, countCoarseNodes = 0;
135 for(
int k = 0; k < this->ghostedNodesPerDir[2]; ++k) {
136 for(
int j = 0; j < this->ghostedNodesPerDir[1]; ++j) {
137 for(
int i = 0; i < this->ghostedNodesPerDir[0]; ++i) {
138 currentIndex = k*this->numGhostedNodes10 + j*this->ghostedNodesPerDir[0] + i;
139 ghostedCoarseNodeCoarseIndices[0] = this->startGhostedCoarseNode[0] + i;
140 ghostedCoarseNodeFineIndices[0] = ghostedCoarseNodeCoarseIndices[0]*this->coarseRate[0];
141 if(ghostedCoarseNodeFineIndices[0] > this->gFineNodesPerDir[0] - 1) {
142 ghostedCoarseNodeFineIndices[0] = this->gFineNodesPerDir[0] - 1;
144 ghostedCoarseNodeCoarseIndices[1] = this->startGhostedCoarseNode[1] + j;
145 ghostedCoarseNodeFineIndices[1] = ghostedCoarseNodeCoarseIndices[1]*this->coarseRate[1];
146 if(ghostedCoarseNodeFineIndices[1] > this->gFineNodesPerDir[1] - 1) {
147 ghostedCoarseNodeFineIndices[1] = this->gFineNodesPerDir[1] - 1;
149 ghostedCoarseNodeCoarseIndices[2] = this->startGhostedCoarseNode[2] + k;
150 ghostedCoarseNodeFineIndices[2] = ghostedCoarseNodeCoarseIndices[2]*this->coarseRate[2];
151 if(ghostedCoarseNodeFineIndices[2] > this->gFineNodesPerDir[2] - 1) {
152 ghostedCoarseNodeFineIndices[2] = this->gFineNodesPerDir[2] - 1;
155 GO myGID = -1, myCoarseGID = -1;
156 LO myLID = -1, myPID = -1, myCoarseLID = -1;
157 getGIDLocalLexicographic(i, j, k, ghostedCoarseNodeFineIndices, myGID, myPID, myLID);
159 int rankIndex = rankIndices[myPID];
160 for(
int dim = 0; dim < 3; ++dim) {
161 if(dim < this->numDimensions) {
162 lCoarseNodeCoarseIndices[dim] = ghostedCoarseNodeCoarseIndices[dim]
163 - coarseMeshData[rankIndex][3 + 2*dim];
166 LO myRankIndexCoarseNodesInDir0 = coarseMeshData[rankIndex][4]
167 - coarseMeshData[rankIndex][3] + 1;
168 LO myRankIndexCoarseNodes10 = (coarseMeshData[rankIndex][6]
169 - coarseMeshData[rankIndex][5] + 1)
170 *myRankIndexCoarseNodesInDir0;
171 myCoarseLID = lCoarseNodeCoarseIndices[2]*myRankIndexCoarseNodes10
172 + lCoarseNodeCoarseIndices[1]*myRankIndexCoarseNodesInDir0
173 + lCoarseNodeCoarseIndices[0];
174 myCoarseGID = myCoarseLID + coarseMeshData[rankIndex][9];
176 ghostedNodeCoarseLIDs[currentIndex] = myCoarseLID;
177 ghostedNodeCoarsePIDs[currentIndex] = myPID;
178 ghostedNodeCoarseGIDs[currentIndex] = myCoarseGID;
180 if(myPID == myRank) {
181 lCoarseNodeCoarseGIDs[countCoarseNodes] = myCoarseGID;
189 template<
class LocalOrdinal,
class GlobalOrdinal,
class Node>
192 Array<GO>& coarseNodeCoarseGIDs,
193 Array<GO>& coarseNodeFineGIDs)
const {
196 coarseNodeCoarseGIDs.resize(this->getNumLocalCoarseNodes());
197 coarseNodeFineGIDs.resize(this->getNumLocalCoarseNodes());
200 ArrayView<const GO> fineNodeGIDs = fineCoordinatesMap->getLocalElementList();
202 Array<GO> coarseStartIndices(3);
203 for(
int dim = 0; dim < 3; ++dim) {
204 coarseStartIndices[dim] = this->coarseMeshData[myRankIndex][2*dim + 3];
209 for(LO coarseLID = 0; coarseLID < this->getNumLocalCoarseNodes(); ++coarseLID) {
210 Array<LO> coarseIndices(3), fineIndices(3), gCoarseIndices(3);
211 this->getCoarseNodeLocalTuple(coarseLID,
215 getCoarseNodeFineLID(coarseIndices[0],coarseIndices[1],coarseIndices[2],fineLID);
216 coarseNodeFineGIDs[coarseLID] = fineNodeGIDs[fineLID];
218 LO myRankIndexCoarseNodesInDir0 = coarseMeshData[myRankIndex][4]
219 - coarseMeshData[myRankIndex][3] + 1;
220 LO myRankIndexCoarseNodes10 = (coarseMeshData[myRankIndex][6]
221 - coarseMeshData[myRankIndex][5] + 1)
222 *myRankIndexCoarseNodesInDir0;
223 LO myCoarseLID = coarseIndices[2]*myRankIndexCoarseNodes10
224 + coarseIndices[1]*myRankIndexCoarseNodesInDir0
226 GO myCoarseGID = myCoarseLID + coarseMeshData[myRankIndex][9];
227 coarseNodeCoarseGIDs[coarseLID] = myCoarseGID;
232 template<
class LocalOrdinal,
class GlobalOrdinal,
class Node>
235 const Array<LO> coarseNodeFineIndices,
236 GO& myGID, LO& myPID, LO& myLID)
const {
238 LO ni = -1, nj = -1, li = -1, lj = -1, lk = -1;
239 LO myRankGuess = myRankIndex;
241 if(iGhosted == 0 && this->ghostInterface[0]) {
243 }
else if((iGhosted == this->ghostedNodesPerDir[0] - 1) && this->ghostInterface[1]) {
246 if(jGhosted == 0 && this->ghostInterface[2]) {
248 }
else if((jGhosted == this->ghostedNodesPerDir[1] - 1) && this->ghostInterface[3]) {
251 if(kGhosted == 0 && this->ghostInterface[4]) {
252 myRankGuess -= pj*pi;
253 }
else if((kGhosted == this->ghostedNodesPerDir[2] - 1) && this->ghostInterface[5]) {
254 myRankGuess += pj*pi;
256 if(coarseNodeFineIndices[0] >= meshData[myRankGuess][3]
257 && coarseNodeFineIndices[0] <= meshData[myRankGuess][4]
258 && coarseNodeFineIndices[1] >= meshData[myRankGuess][5]
259 && coarseNodeFineIndices[1] <= meshData[myRankGuess][6]
260 && coarseNodeFineIndices[2] >= meshData[myRankGuess][7]
261 && coarseNodeFineIndices[2] <= meshData[myRankGuess][8]
262 && myRankGuess < numRanks - 1) {
263 myPID = meshData[myRankGuess][0];
264 ni = meshData[myRankGuess][4] - meshData[myRankGuess][3] + 1;
265 nj = meshData[myRankGuess][6] - meshData[myRankGuess][5] + 1;
266 li = coarseNodeFineIndices[0] - meshData[myRankGuess][3];
267 lj = coarseNodeFineIndices[1] - meshData[myRankGuess][5];
268 lk = coarseNodeFineIndices[2] - meshData[myRankGuess][7];
269 myLID = lk*nj*ni + lj*ni + li;
270 myGID = meshData[myRankGuess][9] + myLID;
274 auto nodeRank = std::find_if(myBlockStart, myBlockEnd,
275 [coarseNodeFineIndices](
const std::vector<GO>& vec){
276 if(coarseNodeFineIndices[0] >= vec[3]
277 && coarseNodeFineIndices[0] <= vec[4]
278 && coarseNodeFineIndices[1] >= vec[5]
279 && coarseNodeFineIndices[1] <= vec[6]
280 && coarseNodeFineIndices[2] >= vec[7]
281 && coarseNodeFineIndices[2] <= vec[8]) {
287 myPID = (*nodeRank)[0];
288 ni = (*nodeRank)[4] - (*nodeRank)[3] + 1;
289 nj = (*nodeRank)[6] - (*nodeRank)[5] + 1;
290 li = coarseNodeFineIndices[0] - (*nodeRank)[3];
291 lj = coarseNodeFineIndices[1] - (*nodeRank)[5];
292 lk = coarseNodeFineIndices[2] - (*nodeRank)[7];
293 myLID = lk*nj*ni + lj*ni + li;
294 myGID = (*nodeRank)[9] + myLID;
298 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
302 std::sort(meshData.begin(), meshData.end(),
303 [](
const std::vector<GO>& a,
const std::vector<GO>& b)->bool {
307 }
else if(a[2] == b[2]) {
310 }
else if(a[7] == b[7]) {
313 }
else if(a[5] == b[5]) {
314 if(a[3] < b[3]) {return true;}
321 numBlocks = meshData[numRanks - 1][2] + 1;
323 myBlockStart = std::lower_bound(meshData.begin(), meshData.end(), myBlock - 1,
324 [] (
const std::vector<GO>& vec,
const GO val)->bool {
325 return (vec[2] < val) ? true : false;
327 myBlockEnd = std::upper_bound(meshData.begin(), meshData.end(), myBlock,
328 [] (
const GO val,
const std::vector<GO>& vec)->bool {
329 return (val < vec[2]) ? true : false;
334 auto myKEnd = std::upper_bound(myBlockStart, myBlockEnd, (*myBlockStart)[3],
335 [] (
const GO val,
const std::vector<GO>& vec)->
bool {
336 return (val < vec[7]) ?
true :
false;
338 auto myJEnd = std::upper_bound(myBlockStart, myKEnd, (*myBlockStart)[3],
339 [] (
const GO val,
const std::vector<GO>& vec)->
bool {
340 return (val < vec[5]) ?
true :
false;
342 pi = std::distance(myBlockStart, myJEnd);
343 pj = std::distance(myBlockStart, myKEnd) / pi;
344 pk = std::distance(myBlockStart, myBlockEnd) / (pj*pi);
347 const int MyRank = myRank;
348 myRankIndex = std::distance(meshData.begin(),
349 std::find_if(myBlockStart, myBlockEnd,
350 [MyRank] (
const std::vector<GO>& vec)->bool {
351 return (vec[0] == MyRank) ? true : false;
356 for(
int rankIndex = 0; rankIndex < numRanks; ++rankIndex) {
357 rankIndices[meshData[rankIndex][0]] = rankIndex;
361 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
364 Array<LO> rankOffset(3);
365 for(
int rank = 0; rank < numRanks; ++rank) {
366 coarseMeshData[rank].resize(10);
367 coarseMeshData[rank][0] = meshData[rank][0];
368 coarseMeshData[rank][1] = meshData[rank][1];
369 coarseMeshData[rank][2] = meshData[rank][2];
370 for(
int dim = 0; dim < 3; ++dim) {
371 coarseMeshData[rank][3 + 2*dim] = meshData[rank][3 + 2*dim] / this->coarseRate[dim];
372 if(meshData[rank][3 + 2*dim] % this->coarseRate[dim] > 0) {
373 ++coarseMeshData[rank][3 + 2*dim];
375 coarseMeshData[rank][3 + 2*dim + 1] = meshData[rank][3 + 2*dim + 1] / this->coarseRate[dim];
376 if(meshData[rank][3 + 2*dim + 1] == this->gFineNodesPerDir[dim] - 1 &&
377 meshData[rank][3 + 2*dim + 1] % this->coarseRate[dim] > 0) {
379 ++coarseMeshData[rank][3 + 2*dim + 1];
383 coarseMeshData[rank][9] = coarseMeshData[rank - 1][9]
384 + (coarseMeshData[rank - 1][8] - coarseMeshData[rank - 1][7] + 1)
385 * (coarseMeshData[rank - 1][6] - coarseMeshData[rank - 1][5] + 1)
386 * (coarseMeshData[rank - 1][4] - coarseMeshData[rank - 1][3] + 1);
391 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
395 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
400 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
404 k = myLID / this->lNumFineNodes10;
405 tmp = myLID % this->lNumFineNodes10;
406 j = tmp / this->lFineNodesPerDir[0];
407 i = tmp % this->lFineNodesPerDir[0];
410 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
414 k = myLID / this->lNumFineNodes10;
415 tmp = myLID % this->lNumFineNodes10;
416 j = tmp / this->lFineNodesPerDir[0];
417 i = tmp % this->lFineNodesPerDir[0];
419 k += this->offsets[2];
420 j += this->offsets[1];
421 i += this->offsets[0];
424 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
429 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
434 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
439 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
443 k = myLID / this->lNumCoarseNodes10;
444 tmp = myLID % this->lNumCoarseNodes10;
445 j = tmp / this->lCoarseNodesPerDir[0];
446 i = tmp % this->lCoarseNodesPerDir[0];
449 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
454 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
459 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
462 myLID = k*this->numGhostedNodes10 + j*this->ghostedNodesPerDir[0] + i;
465 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
470 const LO multiplier[3] = {1, this->lFineNodesPerDir[0], this->lNumFineNodes10};
471 const LO indices[3] = {i, j, k};
474 for(
int dim = 0; dim < 3; ++dim) {
475 if((indices[dim] == this->getLocalCoarseNodesInDir(dim) - 1) && this->meshEdge[2*dim + 1]) {
478 myLID += (this->getLocalFineNodesInDir(dim) - 1)*multiplier[dim];
480 myLID += (indices[dim]*this->getCoarseningRate(dim) + this->getCoarseNodeOffset(dim))
486 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
491 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
Container class for mesh layout and indices calculation.
const bool coupled_
Flag for coupled vs uncoupled aggregation mode, if true aggregation is coupled.
void computeMeshParameters()
Array< GO > startIndices
lowest global tuple (i,j,k) of a node on the local process
Array< int > coarseRate
coarsening rate in each direction
const int numDimensions
Number of spacial dimensions in the problem.
void getGhostedNodesData(const RCP< const Map > fineMap, Array< LO > &ghostedNodeCoarseLIDs, Array< int > &ghostedNodeCoarsePIDs, Array< GO > &ghostedNodeCoarseGIDs) const
const int numRanks
Number of ranks used to decompose the problem.
void computeGlobalCoarseParameters()
int myRankIndex
local process index for record in meshData after sorting.
void sortLocalLexicographicData()
int myBlock
local mesh block ID.
void getGIDLocalLexicographic(const LO iGhosted, const LO jGhosted, const LO kGhosted, const Array< LO > coarseNodeFineIndices, GO &myGID, LO &myPID, LO &myLID) const
LocalLexicographicIndexManager()=default
std::vector< std::vector< GO > > meshData
layout of indices accross all processes.
std::vector< std::vector< GO > > coarseMeshData
layout of indices accross all processes after coarsening.
Array< int > rankIndices
mapping between rank ID and reordered rank ID.
const int myRank
Local rank ID.
void getCoarseNodesData(const RCP< const Map > fineCoordinatesMap, Array< GO > &coarseNodeCoarseGIDs, Array< GO > &coarseNodeFineGIDs) const
void computeCoarseLocalLexicographicData()
Namespace for MueLu classes and methods.