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,
61 :
IndexManager(comm, coupled, false, NumDimensions, interpolationOrder, GFineNodesPerDir, LFineNodesPerDir)
63 , numRanks(NumRanks) {
70 for (
int dim = 0; dim < 3; ++dim) {
72 if (CoarseRate.
size() == 1) {
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>
119 Array<GO>& ghostedNodeCoarseGIDs)
const {
121 ghostedNodeCoarseLIDs.
resize(this->getNumLocalGhostedNodes());
122 ghostedNodeCoarsePIDs.
resize(this->getNumLocalGhostedNodes());
123 ghostedNodeCoarseGIDs.
resize(this->numGhostedNodes);
130 Array<LO> ghostedCoarseNodeCoarseIndices(3), ghostedCoarseNodeFineIndices(3);
132 Array<GO> lCoarseNodeCoarseGIDs(this->lNumCoarseNodes);
133 LO currentIndex = -1, countCoarseNodes = 0;
134 for (
int k = 0; k < this->ghostedNodesPerDir[2]; ++k) {
135 for (
int j = 0; j < this->ghostedNodesPerDir[1]; ++j) {
136 for (
int i = 0; i < this->ghostedNodesPerDir[0]; ++i) {
137 currentIndex = k * this->numGhostedNodes10 + j * this->ghostedNodesPerDir[0] + i;
138 ghostedCoarseNodeCoarseIndices[0] = this->startGhostedCoarseNode[0] + i;
139 ghostedCoarseNodeFineIndices[0] = ghostedCoarseNodeCoarseIndices[0] * this->coarseRate[0];
140 if (ghostedCoarseNodeFineIndices[0] > this->gFineNodesPerDir[0] - 1) {
141 ghostedCoarseNodeFineIndices[0] = this->gFineNodesPerDir[0] - 1;
143 ghostedCoarseNodeCoarseIndices[1] = this->startGhostedCoarseNode[1] + j;
144 ghostedCoarseNodeFineIndices[1] = ghostedCoarseNodeCoarseIndices[1] * this->coarseRate[1];
145 if (ghostedCoarseNodeFineIndices[1] > this->gFineNodesPerDir[1] - 1) {
146 ghostedCoarseNodeFineIndices[1] = this->gFineNodesPerDir[1] - 1;
148 ghostedCoarseNodeCoarseIndices[2] = this->startGhostedCoarseNode[2] + k;
149 ghostedCoarseNodeFineIndices[2] = ghostedCoarseNodeCoarseIndices[2] * this->coarseRate[2];
150 if (ghostedCoarseNodeFineIndices[2] > this->gFineNodesPerDir[2] - 1) {
151 ghostedCoarseNodeFineIndices[2] = this->gFineNodesPerDir[2] - 1;
154 GO myGID = -1, myCoarseGID = -1;
155 LO myLID = -1, myPID = -1, myCoarseLID = -1;
156 getGIDLocalLexicographic(i, j, k, ghostedCoarseNodeFineIndices, myGID, myPID, myLID);
158 int rankIndex = rankIndices[myPID];
159 for (
int dim = 0; dim < 3; ++dim) {
160 if (dim < this->numDimensions) {
161 lCoarseNodeCoarseIndices[dim] = ghostedCoarseNodeCoarseIndices[dim] - coarseMeshData[rankIndex][3 + 2 * dim];
164 LO myRankIndexCoarseNodesInDir0 = coarseMeshData[rankIndex][4] - coarseMeshData[rankIndex][3] + 1;
165 LO myRankIndexCoarseNodes10 = (coarseMeshData[rankIndex][6] - coarseMeshData[rankIndex][5] + 1) * myRankIndexCoarseNodesInDir0;
166 myCoarseLID = lCoarseNodeCoarseIndices[2] * myRankIndexCoarseNodes10 + lCoarseNodeCoarseIndices[1] * myRankIndexCoarseNodesInDir0 + lCoarseNodeCoarseIndices[0];
167 myCoarseGID = myCoarseLID + coarseMeshData[rankIndex][9];
169 ghostedNodeCoarseLIDs[currentIndex] = myCoarseLID;
170 ghostedNodeCoarsePIDs[currentIndex] = myPID;
171 ghostedNodeCoarseGIDs[currentIndex] = myCoarseGID;
173 if (myPID == myRank) {
174 lCoarseNodeCoarseGIDs[countCoarseNodes] = myCoarseGID;
182 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
188 coarseNodeCoarseGIDs.
resize(this->getNumLocalCoarseNodes());
189 coarseNodeFineGIDs.
resize(this->getNumLocalCoarseNodes());
195 for (
int dim = 0; dim < 3; ++dim) {
196 coarseStartIndices[dim] = this->coarseMeshData[myRankIndex][2 * dim + 3];
201 for (
LO coarseLID = 0; coarseLID < this->getNumLocalCoarseNodes(); ++coarseLID) {
202 Array<LO> coarseIndices(3), fineIndices(3), gCoarseIndices(3);
203 this->getCoarseNodeLocalTuple(coarseLID,
207 getCoarseNodeFineLID(coarseIndices[0], coarseIndices[1], coarseIndices[2], fineLID);
208 coarseNodeFineGIDs[coarseLID] = fineNodeGIDs[fineLID];
210 LO myRankIndexCoarseNodesInDir0 = coarseMeshData[myRankIndex][4] - coarseMeshData[myRankIndex][3] + 1;
211 LO myRankIndexCoarseNodes10 = (coarseMeshData[myRankIndex][6] - coarseMeshData[myRankIndex][5] + 1) * myRankIndexCoarseNodesInDir0;
212 LO myCoarseLID = coarseIndices[2] * myRankIndexCoarseNodes10 + coarseIndices[1] * myRankIndexCoarseNodesInDir0 + coarseIndices[0];
213 GO myCoarseGID = myCoarseLID + coarseMeshData[myRankIndex][9];
214 coarseNodeCoarseGIDs[coarseLID] = myCoarseGID;
218 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
222 GO& myGID,
LO& myPID,
LO& myLID)
const {
223 LO ni = -1, nj = -1, li = -1, lj = -1, lk = -1;
224 LO myRankGuess = myRankIndex;
226 if (iGhosted == 0 && this->ghostInterface[0]) {
228 }
else if ((iGhosted == this->ghostedNodesPerDir[0] - 1) && this->ghostInterface[1]) {
231 if (jGhosted == 0 && this->ghostInterface[2]) {
233 }
else if ((jGhosted == this->ghostedNodesPerDir[1] - 1) && this->ghostInterface[3]) {
236 if (kGhosted == 0 && this->ghostInterface[4]) {
237 myRankGuess -= pj * pi;
238 }
else if ((kGhosted == this->ghostedNodesPerDir[2] - 1) && this->ghostInterface[5]) {
239 myRankGuess += pj * pi;
241 if (coarseNodeFineIndices[0] >= meshData[myRankGuess][3] && coarseNodeFineIndices[0] <= meshData[myRankGuess][4] && coarseNodeFineIndices[1] >= meshData[myRankGuess][5] && coarseNodeFineIndices[1] <= meshData[myRankGuess][6] && coarseNodeFineIndices[2] >= meshData[myRankGuess][7] && coarseNodeFineIndices[2] <= meshData[myRankGuess][8] && myRankGuess < numRanks - 1) {
242 myPID = meshData[myRankGuess][0];
243 ni = meshData[myRankGuess][4] - meshData[myRankGuess][3] + 1;
244 nj = meshData[myRankGuess][6] - meshData[myRankGuess][5] + 1;
245 li = coarseNodeFineIndices[0] - meshData[myRankGuess][3];
246 lj = coarseNodeFineIndices[1] - meshData[myRankGuess][5];
247 lk = coarseNodeFineIndices[2] - meshData[myRankGuess][7];
248 myLID = lk * nj * ni + lj * ni + li;
249 myGID = meshData[myRankGuess][9] + myLID;
253 auto nodeRank = std::find_if(myBlockStart, myBlockEnd,
254 [coarseNodeFineIndices](
const std::vector<GO>& vec) {
255 if (coarseNodeFineIndices[0] >= vec[3] && coarseNodeFineIndices[0] <= vec[4] && coarseNodeFineIndices[1] >= vec[5] && coarseNodeFineIndices[1] <= vec[6] && coarseNodeFineIndices[2] >= vec[7] && coarseNodeFineIndices[2] <= vec[8]) {
261 myPID = (*nodeRank)[0];
262 ni = (*nodeRank)[4] - (*nodeRank)[3] + 1;
263 nj = (*nodeRank)[6] - (*nodeRank)[5] + 1;
264 li = coarseNodeFineIndices[0] - (*nodeRank)[3];
265 lj = coarseNodeFineIndices[1] - (*nodeRank)[5];
266 lk = coarseNodeFineIndices[2] - (*nodeRank)[7];
267 myLID = lk * nj * ni + lj * ni + li;
268 myGID = (*nodeRank)[9] + myLID;
272 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
275 std::sort(meshData.begin(), meshData.end(),
276 [](
const std::vector<GO>& a,
const std::vector<GO>& b) ->
bool {
280 }
else if (a[2] == b[2]) {
283 }
else if (a[7] == b[7]) {
286 }
else if (a[5] == b[5]) {
296 numBlocks = meshData[numRanks - 1][2] + 1;
298 myBlockStart = std::lower_bound(meshData.begin(), meshData.end(), myBlock - 1,
299 [](
const std::vector<GO>& vec,
const GO val) ->
bool {
300 return (vec[2] < val) ?
true :
false;
302 myBlockEnd = std::upper_bound(meshData.begin(), meshData.end(), myBlock,
303 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
304 return (val < vec[2]) ?
true :
false;
309 auto myKEnd = std::upper_bound(myBlockStart, myBlockEnd, (*myBlockStart)[3],
310 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
311 return (val < vec[7]) ?
true :
false;
313 auto myJEnd = std::upper_bound(myBlockStart, myKEnd, (*myBlockStart)[3],
314 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
315 return (val < vec[5]) ?
true :
false;
317 pi = std::distance(myBlockStart, myJEnd);
318 pj = std::distance(myBlockStart, myKEnd) / pi;
319 pk = std::distance(myBlockStart, myBlockEnd) / (pj * pi);
322 const int MyRank = myRank;
323 myRankIndex = std::distance(meshData.begin(),
324 std::find_if(myBlockStart, myBlockEnd,
325 [MyRank](
const std::vector<GO>& vec) ->
bool {
326 return (vec[0] == MyRank) ?
true :
false;
330 for (
int rankIndex = 0; rankIndex < numRanks; ++rankIndex) {
331 rankIndices[meshData[rankIndex][0]] = rankIndex;
335 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
339 for (
int rank = 0; rank < numRanks; ++rank) {
340 coarseMeshData[rank].resize(10);
341 coarseMeshData[rank][0] = meshData[rank][0];
342 coarseMeshData[rank][1] = meshData[rank][1];
343 coarseMeshData[rank][2] = meshData[rank][2];
344 for (
int dim = 0; dim < 3; ++dim) {
345 coarseMeshData[rank][3 + 2 * dim] = meshData[rank][3 + 2 * dim] / this->coarseRate[dim];
346 if (meshData[rank][3 + 2 * dim] % this->coarseRate[dim] > 0) {
347 ++coarseMeshData[rank][3 + 2 * dim];
349 coarseMeshData[rank][3 + 2 * dim + 1] = meshData[rank][3 + 2 * dim + 1] / this->coarseRate[dim];
350 if (meshData[rank][3 + 2 * dim + 1] == this->gFineNodesPerDir[dim] - 1 &&
351 meshData[rank][3 + 2 * dim + 1] % this->coarseRate[dim] > 0) {
353 ++coarseMeshData[rank][3 + 2 * dim + 1];
357 coarseMeshData[rank][9] = coarseMeshData[rank - 1][9] + (coarseMeshData[rank - 1][8] - coarseMeshData[rank - 1][7] + 1) * (coarseMeshData[rank - 1][6] - coarseMeshData[rank - 1][5] + 1) * (coarseMeshData[rank - 1][4] - coarseMeshData[rank - 1][3] + 1);
362 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
366 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
371 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
375 k = myLID / this->lNumFineNodes10;
376 tmp = myLID % this->lNumFineNodes10;
377 j = tmp / this->lFineNodesPerDir[0];
378 i = tmp % this->lFineNodesPerDir[0];
381 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
385 k = myLID / this->lNumFineNodes10;
386 tmp = myLID % this->lNumFineNodes10;
387 j = tmp / this->lFineNodesPerDir[0];
388 i = tmp % this->lFineNodesPerDir[0];
390 k += this->offsets[2];
391 j += this->offsets[1];
392 i += this->offsets[0];
395 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
400 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
405 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
410 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
414 k = myLID / this->lNumCoarseNodes10;
415 tmp = myLID % this->lNumCoarseNodes10;
416 j = tmp / this->lCoarseNodesPerDir[0];
417 i = tmp % this->lCoarseNodesPerDir[0];
420 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
425 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
430 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
433 myLID = k * this->numGhostedNodes10 + j * this->ghostedNodesPerDir[0] + i;
436 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
441 const LO multiplier[3] = {1, this->lFineNodesPerDir[0], this->lNumFineNodes10};
442 const LO indices[3] = {i, j, k};
445 for (
int dim = 0; dim < 3; ++dim) {
446 if ((indices[dim] == this->getLocalCoarseNodesInDir(dim) - 1) && this->meshEdge[2 * dim + 1]) {
449 myLID += (this->getLocalFineNodesInDir(dim) - 1) * multiplier[dim];
451 myLID += (indices[dim] * this->getCoarseningRate(dim) + this->getCoarseNodeOffset(dim)) * multiplier[dim];
456 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
461 template <
class LocalOrdinal,
class GlobalOrdinal,
class Node>
const bool coupled_
Flag for coupled vs uncoupled aggregation mode, if true aggregation is coupled.
void getCoarseNodeGlobalTuple(const GO myGID, GO &i, GO &j, GO &k) const
void getCoarseNodeLocalTuple(const LO myLID, LO &i, LO &j, LO &k) const
void getGhostedNodesData(const RCP< const Map > fineMap, Array< LO > &ghostedNodeCoarseLIDs, Array< int > &ghostedNodeCoarsePIDs, Array< GO > &ghostedNodeCoarseGIDs) const
void getFineNodeLocalTuple(const LO myLID, LO &i, LO &j, LO &k) const
std::vector< std::vector< GO > > getCoarseMeshData() const
void sortLocalLexicographicData()
void getFineNodeGhostedTuple(const LO myLID, LO &i, LO &j, LO &k) const
void getGhostedNodeCoarseLID(const LO i, const LO j, const LO k, LO &myLID) const
void getGhostedNodeFineLID(const LO i, const LO j, const LO k, LO &myLID) const
int myRankIndex
local process index for record in meshData after sorting.
void getCoarseNodesData(const RCP< const Map > fineCoordinatesMap, Array< GO > &coarseNodeCoarseGIDs, Array< GO > &coarseNodeFineGIDs) const
Array< GO > startIndices
lowest global tuple (i,j,k) of a node on the local process
void getCoarseNodeFineLID(const LO i, const LO j, const LO k, LO &myLID) const
std::vector< std::vector< GO > > meshData
layout of indices accross all processes.
void getCoarseNodeLID(const LO i, const LO j, const LO k, LO &myLID) const
int myBlock
local mesh block ID.
Array< int > coarseRate
coarsening rate in each direction
void getFineNodeGlobalTuple(const GO myGID, GO &i, GO &j, GO &k) const
const int numDimensions
Number of spacial dimensions in the problem.
Array< int > rankIndices
mapping between rank ID and reordered rank ID.
void getFineNodeGID(const GO i, const GO j, const GO k, GO &myGID) const
const int numRanks
Number of ranks used to decompose the problem.
void computeGlobalCoarseParameters()
void resize(size_type new_size, const value_type &x=value_type())
void getCoarseNodeGID(const GO i, const GO j, const GO k, GO &myGID) const
void getFineNodeLID(const LO i, const LO j, const LO k, LO &myLID) const
LocalLexicographicIndexManager()=default
void computeCoarseLocalLexicographicData()
void getCoarseNodeGhostedLID(const LO i, const LO j, const LO k, LO &myLID) const
void getGIDLocalLexicographic(const LO iGhosted, const LO jGhosted, const LO kGhosted, const Array< LO > coarseNodeFineIndices, GO &myGID, LO &myPID, LO &myLID) const
void computeMeshParameters()
const int myRank
Local rank ID.
Container class for mesh layout and indices calculation.
std::vector< std::vector< GO > > coarseMeshData
layout of indices accross all processes after coarsening.