46 #ifndef MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
47 #define MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
57 #include <Xpetra_CrsMatrixWrap.hpp>
60 #include <Xpetra_MapFactory.hpp>
61 #include <Xpetra_MultiVectorFactory.hpp>
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
81 validParamList->
set<std::string>(
"Coarsen",
"{2}",
82 "Coarsening rate in all spatial dimensions");
83 validParamList->
set<
int>(
"order", 1,
84 "Order of the interpolation scheme used");
86 "Generating factory of the matrix A");
88 "Generating factory of the nullspace");
90 "Generating factory for coorindates");
92 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
94 "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
95 validParamList->
set<std::string>(
"meshLayout",
"Global Lexicographic",
96 "Type of mesh ordering");
98 "Mesh ordering associated data");
100 return validParamList;
103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 Input(fineLevel,
"A");
107 Input(fineLevel,
"Nullspace");
108 Input(fineLevel,
"Coordinates");
116 "gNodesPerDim was not provided by the user on level0!");
119 Input(fineLevel,
"gNodesPerDim");
129 "lNodesPerDim was not provided by the user on level0!");
132 Input(fineLevel,
"lNodesPerDim");
136 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
138 Level& coarseLevel)
const {
139 return BuildP(fineLevel, coarseLevel);
142 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 Level& coarseLevel)
const {
150 RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel,
"Nullspace");
151 RCP<xdMV> fineCoords = Get<RCP<xdMV> >(fineLevel,
"Coordinates");
158 const LO blkSize = A->GetFixedBlockSize();
163 myGeometry->meshLayout = pL.
get<std::string>(
"meshLayout");
164 if (fineLevel.GetLevelID() == 0) {
165 if (myGeometry->meshLayout ==
"Local Lexicographic") {
169 "The meshData array is empty, somehow the input for geometric"
170 " multigrid are not captured correctly.");
171 myGeometry->meshData.resize(rowMap->getComm()->getSize());
172 for (
int i = 0; i < rowMap->getComm()->getSize(); ++i) {
173 myGeometry->meshData[i].resize(10);
174 for (
int j = 0; j < 10; ++j) {
175 myGeometry->meshData[i][j] = meshData[10 * i + j];
182 "Coordinates cannot be accessed from fine level!");
183 myGeometry->numDimensions = fineCoords->getNumVectors();
186 if (fineLevel.GetLevelID() == 0) {
191 myGeometry->gFineNodesPerDir = Get<Array<GO> >(fineLevel,
"gNodesPerDim");
194 myGeometry->lFineNodesPerDir = Get<Array<LO> >(fineLevel,
"lNodesPerDim");
196 myGeometry->gNumFineNodes10 = myGeometry->gFineNodesPerDir[1] * myGeometry->gFineNodesPerDir[0];
197 myGeometry->gNumFineNodes = myGeometry->gFineNodesPerDir[2] * myGeometry->gNumFineNodes10;
198 myGeometry->lNumFineNodes10 = myGeometry->lFineNodesPerDir[1] * myGeometry->lFineNodesPerDir[0];
199 myGeometry->lNumFineNodes = myGeometry->lFineNodesPerDir[2] * myGeometry->lNumFineNodes10;
203 "The local number of elements in Coordinates is not equal to the"
204 " number of nodes given by: lNodesPerDim!");
207 "The global number of elements in Coordinates is not equal to the"
208 " number of nodes given by: gNodesPerDim!");
211 std::string coarsenRate = pL.
get<std::string>(
"Coarsen");
214 crates = Teuchos::fromStringToArray<LO>(coarsenRate);
216 GetOStream(
Errors, -1) <<
" *** Coarsen must be a string convertible into an array! *** "
222 "Coarsen must have at least as many components as the number of"
223 " spatial dimensions in the problem.");
225 for (LO i = 0; i < 3; ++i) {
226 if (i < myGeometry->numDimensions) {
227 if (crates.size() == 1) {
228 myGeometry->coarseRate[i] = crates[0];
229 }
else if (crates.size() == myGeometry->numDimensions) {
230 myGeometry->coarseRate[i] = crates[i];
233 myGeometry->coarseRate[i] = 1;
237 int interpolationOrder = pL.
get<
int>(
"order");
240 "The interpolation order can only be set to 0 or 1.");
244 for (LO i = 0; i < myGeometry->numDimensions; ++i) {
247 for (LO i = 0; i < myGeometry->numDimensions; ++i) {
248 TEUCHOS_TEST_FOR_EXCEPTION(mapDirG2L[i] > myGeometry->numDimensions,
250 "axis permutation values must all be less than"
251 " the number of spatial dimensions.");
252 mapDirL2G[mapDirG2L[i]] = i;
259 if ((fineLevel.GetLevelID() == 0) && (myGeometry->meshLayout !=
"Global Lexicographic")) {
260 MeshLayoutInterface(interpolationOrder, blkSize, fineCoordsMap, myGeometry,
261 ghostedCoarseNodes, lCoarseNodesGIDs);
266 GetCoarsePoints(interpolationOrder, blkSize, fineCoordsMap, myGeometry, ghostedCoarseNodes,
281 GO lnnzP = Teuchos::as<LO>(std::pow(interpolationOrder + 1, myGeometry->numDimensions)) * (myGeometry->lNumFineNodes - myGeometry->lNumCoarseNodes) + myGeometry->lNumCoarseNodes;
282 lnnzP = lnnzP * blkSize;
283 GO gnnzP = Teuchos::as<LO>(std::pow(interpolationOrder + 1, myGeometry->numDimensions)) * (myGeometry->gNumFineNodes - myGeometry->gNumCoarseNodes) + myGeometry->gNumCoarseNodes;
284 gnnzP = gnnzP * blkSize;
286 GetOStream(
Runtime1) <<
"P size = " << blkSize * myGeometry->gNumFineNodes
287 <<
" x " << blkSize * myGeometry->gNumCoarseNodes << std::endl;
288 GetOStream(
Runtime1) <<
"P Fine grid : " << myGeometry->gFineNodesPerDir[0] <<
" -- "
289 << myGeometry->gFineNodesPerDir[1] <<
" -- "
290 << myGeometry->gFineNodesPerDir[2] << std::endl;
291 GetOStream(
Runtime1) <<
"P Coarse grid : " << myGeometry->gCoarseNodesPerDir[0] <<
" -- "
292 << myGeometry->gCoarseNodesPerDir[1] <<
" -- "
293 << myGeometry->gCoarseNodesPerDir[2] << std::endl;
294 GetOStream(
Runtime1) <<
"P nnz estimate: " << gnnzP << std::endl;
296 MakeGeneralGeometricP(myGeometry, fineCoords, lnnzP, blkSize, stridedDomainMapP,
297 A, P, coarseCoords, ghostedCoarseNodes, lCoarseNodesGIDs,
301 if (A->IsView(
"stridedMaps") ==
true) {
302 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMapP);
304 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMapP);
308 Set(coarseLevel,
"P", P);
309 Set(coarseLevel,
"coarseCoordinates", coarseCoords);
310 Set<Array<GO> >(coarseLevel,
"gCoarseNodesPerDim", myGeometry->gCoarseNodesPerDir);
311 Set<Array<LO> >(coarseLevel,
"lCoarseNodesPerDim", myGeometry->lCoarseNodesPerDir);
315 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
316 fineNullspace->getNumVectors());
319 Set(coarseLevel,
"Nullspace", coarseNullspace);
322 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
333 int numRanks = fineCoordsMap->getComm()->getSize();
334 int myRank = fineCoordsMap->getComm()->getRank();
336 myGeo->myBlock = myGeo->meshData[myRank][2];
337 myGeo->startIndices[0] = myGeo->meshData[myRank][3];
338 myGeo->startIndices[1] = myGeo->meshData[myRank][5];
339 myGeo->startIndices[2] = myGeo->meshData[myRank][7];
340 myGeo->startIndices[3] = myGeo->meshData[myRank][4];
341 myGeo->startIndices[4] = myGeo->meshData[myRank][6];
342 myGeo->startIndices[5] = myGeo->meshData[myRank][8];
343 std::sort(myGeo->meshData.begin(), myGeo->meshData.end(),
344 [](
const std::vector<GO>& a,
const std::vector<GO>& b) ->
bool {
348 }
else if (a[2] == b[2]) {
351 }
else if (a[7] == b[7]) {
354 }
else if (a[5] == b[5]) {
364 myGeo->numBlocks = myGeo->meshData[numRanks - 1][2] + 1;
366 auto myBlockStart = std::lower_bound(myGeo->meshData.begin(), myGeo->meshData.end(),
368 [](
const std::vector<GO>& vec,
const GO val) ->
bool {
369 return (vec[2] < val) ?
true :
false;
371 auto myBlockEnd = std::upper_bound(myGeo->meshData.begin(), myGeo->meshData.end(),
373 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
374 return (val < vec[2]) ?
true :
false;
379 auto myKEnd = std::upper_bound(myBlockStart, myBlockEnd, (*myBlockStart)[3],
380 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
381 return (val < vec[7]) ?
true :
false;
383 auto myJEnd = std::upper_bound(myBlockStart, myKEnd, (*myBlockStart)[3],
384 [](
const GO val,
const std::vector<GO>& vec) ->
bool {
385 return (val < vec[5]) ?
true :
false;
387 LO pi = std::distance(myBlockStart, myJEnd);
388 LO pj = std::distance(myBlockStart, myKEnd) / pi;
389 LO pk = std::distance(myBlockStart, myBlockEnd) / (pj * pi);
392 LO myRankIndex = std::distance(myGeo->meshData.begin(),
393 std::find_if(myBlockStart, myBlockEnd,
394 [myRank](
const std::vector<GO>& vec) ->
bool {
395 return (vec[0] == myRank) ?
true :
false;
398 for (
int dim = 0; dim < 3; ++dim) {
399 if (dim < myGeo->numDimensions) {
400 myGeo->offsets[dim] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
401 myGeo->offsets[dim + 3] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
407 for (
int dim = 0; dim < 3; ++dim) {
408 if (dim < myGeo->numDimensions && (myGeo->startIndices[dim] % myGeo->coarseRate[dim] != 0 ||
409 myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1)) {
410 myGeo->ghostInterface[2 * dim] =
true;
412 if (dim < myGeo->numDimensions && myGeo->startIndices[dim + 3] != myGeo->gFineNodesPerDir[dim] - 1 && (myGeo->lFineNodesPerDir[dim] == 1 || myGeo->startIndices[dim + 3] % myGeo->coarseRate[dim] != 0)) {
413 myGeo->ghostInterface[2 * dim + 1] =
true;
429 for (
int i = 0; i < 3; ++i) {
430 if (i < myGeo->numDimensions) {
433 myGeo->gCoarseNodesPerDir[i] = (myGeo->gFineNodesPerDir[i] - 1) / myGeo->coarseRate[i];
434 myGeo->endRate[i] = Teuchos::as<LO>((myGeo->gFineNodesPerDir[i] - 1) % myGeo->coarseRate[i]);
435 if (myGeo->endRate[i] == 0) {
436 myGeo->endRate[i] = myGeo->coarseRate[i];
437 ++myGeo->gCoarseNodesPerDir[i];
439 myGeo->gCoarseNodesPerDir[i] += 2;
442 myGeo->endRate[i] = 1;
443 myGeo->gCoarseNodesPerDir[i] = 1;
447 myGeo->gNumCoarseNodes = myGeo->gCoarseNodesPerDir[0] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[2];
449 for (
LO i = 0; i < 3; ++i) {
450 if (i < myGeo->numDimensions) {
454 if ((myGeo->startIndices[i] + myGeo->lFineNodesPerDir[i]) == myGeo->gFineNodesPerDir[i]) {
455 myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] - myGeo->endRate[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i] + 1;
456 if (myGeo->offsets[i] == 0) {
457 ++myGeo->lCoarseNodesPerDir[i];
460 myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i];
461 if (myGeo->offsets[i] == 0) {
462 ++myGeo->lCoarseNodesPerDir[i];
466 myGeo->lCoarseNodesPerDir[i] = 1;
470 if (myGeo->lFineNodesPerDir[i] < 1) {
471 myGeo->lCoarseNodesPerDir[i] = 0;
478 myGeo->lNumCoarseNodes = myGeo->lCoarseNodesPerDir[0] * myGeo->lCoarseNodesPerDir[1] * myGeo->lCoarseNodesPerDir[2];
481 for (
int dim = 0; dim < 3; ++dim) {
485 if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1 &&
486 myGeo->startIndices[dim] % myGeo->coarseRate[dim] == 0) {
487 myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim] - 1;
489 myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim];
491 myGeo->ghostedCoarseNodesPerDir[dim] = myGeo->lCoarseNodesPerDir[dim];
493 if (myGeo->ghostInterface[2 * dim]) {
494 myGeo->ghostedCoarseNodesPerDir[dim] += 1;
497 if (myGeo->ghostInterface[2 * dim + 1]) {
498 myGeo->ghostedCoarseNodesPerDir[dim] += 1;
501 myGeo->lNumGhostedNodes = myGeo->ghostedCoarseNodesPerDir[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0];
502 myGeo->lNumGhostNodes = myGeo->lNumGhostedNodes - myGeo->lNumCoarseNodes;
503 ghostedCoarseNodes->PIDs.resize(myGeo->lNumGhostedNodes);
504 ghostedCoarseNodes->LIDs.resize(myGeo->lNumGhostedNodes);
505 ghostedCoarseNodes->GIDs.resize(myGeo->lNumGhostedNodes);
506 ghostedCoarseNodes->coarseGIDs.resize(myGeo->lNumGhostedNodes);
507 ghostedCoarseNodes->colInds.resize(myGeo->lNumGhostedNodes);
508 lCoarseNodesGIDs[0].resize(myGeo->lNumCoarseNodes);
509 lCoarseNodesGIDs[1].resize(myGeo->lNumCoarseNodes);
516 Array<LO> coarseNodeCoarseIndices(3), coarseNodeFineIndices(3);
517 LO currentIndex = -1, countCoarseNodes = 0;
518 for (
int k = 0; k < myGeo->ghostedCoarseNodesPerDir[2]; ++k) {
519 for (
int j = 0; j < myGeo->ghostedCoarseNodesPerDir[1]; ++j) {
520 for (
int i = 0; i < myGeo->ghostedCoarseNodesPerDir[0]; ++i) {
521 currentIndex = k * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + j * myGeo->ghostedCoarseNodesPerDir[0] + i;
522 coarseNodeCoarseIndices[0] = myGeo->startGhostedCoarseNode[0] + i;
523 coarseNodeFineIndices[0] = coarseNodeCoarseIndices[0] * myGeo->coarseRate[0];
524 if (coarseNodeFineIndices[0] > myGeo->gFineNodesPerDir[0] - 1) {
525 coarseNodeFineIndices[0] = myGeo->gFineNodesPerDir[0] - 1;
527 coarseNodeCoarseIndices[1] = myGeo->startGhostedCoarseNode[1] + j;
528 coarseNodeFineIndices[1] = coarseNodeCoarseIndices[1] * myGeo->coarseRate[1];
529 if (coarseNodeFineIndices[1] > myGeo->gFineNodesPerDir[1] - 1) {
530 coarseNodeFineIndices[1] = myGeo->gFineNodesPerDir[1] - 1;
532 coarseNodeCoarseIndices[2] = myGeo->startGhostedCoarseNode[2] + k;
533 coarseNodeFineIndices[2] = coarseNodeCoarseIndices[2] * myGeo->coarseRate[2];
534 if (coarseNodeFineIndices[2] > myGeo->gFineNodesPerDir[2] - 1) {
535 coarseNodeFineIndices[2] = myGeo->gFineNodesPerDir[2] - 1;
537 GO myGID = -1, myCoarseGID = -1;
538 LO myLID = -1, myPID = -1;
539 GetGIDLocalLexicographic(i, j, k, coarseNodeFineIndices, myGeo, myRankIndex, pi, pj, pk,
540 myBlockStart, myBlockEnd, myGID, myPID, myLID);
541 myCoarseGID = coarseNodeCoarseIndices[0] + coarseNodeCoarseIndices[1] * myGeo->gCoarseNodesPerDir[0] + coarseNodeCoarseIndices[2] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[0];
542 ghostedCoarseNodes->PIDs[currentIndex] = myPID;
543 ghostedCoarseNodes->LIDs[currentIndex] = myLID;
544 ghostedCoarseNodes->GIDs[currentIndex] = myGID;
545 ghostedCoarseNodes->coarseGIDs[currentIndex] = myCoarseGID;
546 if (myPID == myRank) {
547 lCoarseNodesGIDs[0][countCoarseNodes] = myCoarseGID;
548 lCoarseNodesGIDs[1][countCoarseNodes] = myGID;
557 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
570 myGeo->startIndices[2] = fineCoordsMap->getMinGlobalIndex() / (myGeo->gFineNodesPerDir[1] * myGeo->gFineNodesPerDir[0]);
571 tmp = fineCoordsMap->getMinGlobalIndex() % (myGeo->gFineNodesPerDir[1] * myGeo->gFineNodesPerDir[0]);
572 myGeo->startIndices[1] = tmp / myGeo->gFineNodesPerDir[0];
573 myGeo->startIndices[0] = tmp % myGeo->gFineNodesPerDir[0];
575 for (
int dim = 0; dim < 3; ++dim) {
576 myGeo->startIndices[dim + 3] = myGeo->startIndices[dim] + myGeo->lFineNodesPerDir[dim] - 1;
579 for (
int dim = 0; dim < 3; ++dim) {
580 if (dim < myGeo->numDimensions) {
581 myGeo->offsets[dim] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
582 myGeo->offsets[dim + 3] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
588 for (
int dim = 0; dim < 3; ++dim) {
589 if (dim < myGeo->numDimensions && (myGeo->startIndices[dim] % myGeo->coarseRate[dim] != 0 ||
590 myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1)) {
591 myGeo->ghostInterface[2 * dim] =
true;
593 if (dim < myGeo->numDimensions && myGeo->startIndices[dim + 3] != myGeo->gFineNodesPerDir[dim] - 1 && (myGeo->lFineNodesPerDir[dim] == 1 || myGeo->startIndices[dim + 3] % myGeo->coarseRate[dim] != 0)) {
594 myGeo->ghostInterface[2 * dim + 1] =
true;
610 for (
int i = 0; i < 3; ++i) {
611 if (i < myGeo->numDimensions) {
614 myGeo->gCoarseNodesPerDir[i] = (myGeo->gFineNodesPerDir[i] - 1) / myGeo->coarseRate[i];
615 myGeo->endRate[i] = Teuchos::as<LO>((myGeo->gFineNodesPerDir[i] - 1) % myGeo->coarseRate[i]);
616 if (myGeo->endRate[i] == 0) {
617 myGeo->endRate[i] = myGeo->coarseRate[i];
618 ++myGeo->gCoarseNodesPerDir[i];
620 myGeo->gCoarseNodesPerDir[i] += 2;
623 myGeo->endRate[i] = 1;
624 myGeo->gCoarseNodesPerDir[i] = 1;
628 myGeo->gNumCoarseNodes = myGeo->gCoarseNodesPerDir[0] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[2];
630 for (
LO i = 0; i < 3; ++i) {
631 if (i < myGeo->numDimensions) {
635 if ((myGeo->startIndices[i] + myGeo->lFineNodesPerDir[i]) == myGeo->gFineNodesPerDir[i]) {
636 myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] - myGeo->endRate[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i] + 1;
637 if (myGeo->offsets[i] == 0) {
638 ++myGeo->lCoarseNodesPerDir[i];
641 myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i];
642 if (myGeo->offsets[i] == 0) {
643 ++myGeo->lCoarseNodesPerDir[i];
647 myGeo->lCoarseNodesPerDir[i] = 1;
651 if (myGeo->lFineNodesPerDir[i] < 1) {
652 myGeo->lCoarseNodesPerDir[i] = 0;
659 myGeo->lNumCoarseNodes = myGeo->lCoarseNodesPerDir[0] * myGeo->lCoarseNodesPerDir[1] * myGeo->lCoarseNodesPerDir[2];
662 bool ghostedDir[6] = {
false};
663 for (
int dim = 0; dim < 3; ++dim) {
667 if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1 &&
668 myGeo->startIndices[dim] % myGeo->coarseRate[dim] == 0) {
669 myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim] - 1;
671 myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim];
673 myGeo->ghostedCoarseNodesPerDir[dim] = myGeo->lCoarseNodesPerDir[dim];
675 if (myGeo->ghostInterface[2 * dim]) {
676 myGeo->ghostedCoarseNodesPerDir[dim] += 1;
677 ghostedDir[2 * dim] =
true;
680 if (myGeo->ghostInterface[2 * dim + 1]) {
681 myGeo->ghostedCoarseNodesPerDir[dim] += 1;
682 ghostedDir[2 * dim + 1] =
true;
685 myGeo->lNumGhostedNodes = myGeo->ghostedCoarseNodesPerDir[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0];
686 myGeo->lNumGhostNodes = myGeo->lNumGhostedNodes - myGeo->lNumCoarseNodes;
687 ghostedCoarseNodes->PIDs.resize(myGeo->lNumGhostedNodes);
688 ghostedCoarseNodes->LIDs.resize(myGeo->lNumGhostedNodes);
689 ghostedCoarseNodes->GIDs.resize(myGeo->lNumGhostedNodes);
690 ghostedCoarseNodes->coarseGIDs.resize(myGeo->lNumGhostedNodes);
691 ghostedCoarseNodes->colInds.resize(myGeo->lNumGhostedNodes);
692 lCoarseNodesGIDs[0].resize(myGeo->lNumCoarseNodes);
693 lCoarseNodesGIDs[1].resize(myGeo->lNumCoarseNodes);
700 Array<LO> coarseNodeCoarseIndices(3), coarseNodeFineIndices(3), ijk(3);
701 LO currentIndex = -1, countCoarseNodes = 0;
702 for (ijk[2] = 0; ijk[2] < myGeo->ghostedCoarseNodesPerDir[2]; ++ijk[2]) {
703 for (ijk[1] = 0; ijk[1] < myGeo->ghostedCoarseNodesPerDir[1]; ++ijk[1]) {
704 for (ijk[0] = 0; ijk[0] < myGeo->ghostedCoarseNodesPerDir[0]; ++ijk[0]) {
705 currentIndex = ijk[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + ijk[1] * myGeo->ghostedCoarseNodesPerDir[0] + ijk[0];
706 coarseNodeCoarseIndices[0] = myGeo->startGhostedCoarseNode[0] + ijk[0];
707 coarseNodeFineIndices[0] = coarseNodeCoarseIndices[0] * myGeo->coarseRate[0];
708 if (coarseNodeFineIndices[0] > myGeo->gFineNodesPerDir[0] - 1) {
709 coarseNodeFineIndices[0] = myGeo->gFineNodesPerDir[0] - 1;
711 coarseNodeCoarseIndices[1] = myGeo->startGhostedCoarseNode[1] + ijk[1];
712 coarseNodeFineIndices[1] = coarseNodeCoarseIndices[1] * myGeo->coarseRate[1];
713 if (coarseNodeFineIndices[1] > myGeo->gFineNodesPerDir[1] - 1) {
714 coarseNodeFineIndices[1] = myGeo->gFineNodesPerDir[1] - 1;
716 coarseNodeCoarseIndices[2] = myGeo->startGhostedCoarseNode[2] + ijk[2];
717 coarseNodeFineIndices[2] = coarseNodeCoarseIndices[2] * myGeo->coarseRate[2];
718 if (coarseNodeFineIndices[2] > myGeo->gFineNodesPerDir[2] - 1) {
719 coarseNodeFineIndices[2] = myGeo->gFineNodesPerDir[2] - 1;
721 GO myGID = 0, myCoarseGID = -1;
723 factor[2] = myGeo->gNumFineNodes10;
724 factor[1] = myGeo->gFineNodesPerDir[0];
726 for (
int dim = 0; dim < 3; ++dim) {
727 if (dim < myGeo->numDimensions) {
728 if (myGeo->startIndices[dim] - myGeo->offsets[dim] + ijk[dim] * myGeo->coarseRate[dim] < myGeo->gFineNodesPerDir[dim] - 1) {
729 myGID += (myGeo->startIndices[dim] - myGeo->offsets[dim] + ijk[dim] * myGeo->coarseRate[dim]) * factor[dim];
731 myGID += (myGeo->startIndices[dim] - myGeo->offsets[dim] + (ijk[dim] - 1) * myGeo->coarseRate[dim] + myGeo->endRate[dim]) * factor[dim];
735 myCoarseGID = coarseNodeCoarseIndices[0] + coarseNodeCoarseIndices[1] * myGeo->gCoarseNodesPerDir[0] + coarseNodeCoarseIndices[2] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[0];
736 ghostedCoarseNodes->GIDs[currentIndex] = myGID;
737 ghostedCoarseNodes->coarseGIDs[currentIndex] = myCoarseGID;
738 if ((!ghostedDir[0] || ijk[0] != 0) && (!ghostedDir[2] || ijk[1] != 0) && (!ghostedDir[4] || ijk[2] != 0) && (!ghostedDir[1] || ijk[0] != myGeo->ghostedCoarseNodesPerDir[0] - 1) && (!ghostedDir[3] || ijk[1] != myGeo->ghostedCoarseNodesPerDir[1] - 1) && (!ghostedDir[5] || ijk[2] != myGeo->ghostedCoarseNodesPerDir[2] - 1)) {
739 lCoarseNodesGIDs[0][countCoarseNodes] = myCoarseGID;
740 lCoarseNodesGIDs[1][countCoarseNodes] = myGID;
746 Array<int> ghostsPIDs(myGeo->lNumGhostedNodes);
747 Array<LO> ghostsLIDs(myGeo->lNumGhostedNodes);
748 fineCoordsMap->getRemoteIndexList(ghostedCoarseNodes->GIDs(),
749 ghostedCoarseNodes->PIDs(),
750 ghostedCoarseNodes->LIDs());
753 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
757 const LO nnzP,
const LO dofsPerNode,
761 int interpolationOrder)
const {
783 LO myRank = Amat->getRowMap()->getComm()->getRank();
784 GO numGloCols = dofsPerNode * myGeo->gNumCoarseNodes;
804 std::vector<NodeID> colMapOrdering(myGeo->lNumGhostedNodes);
805 for (LO ind = 0; ind < myGeo->lNumGhostedNodes; ++ind) {
806 colMapOrdering[ind].GID = ghostedCoarseNodes->GIDs[ind];
807 if (ghostedCoarseNodes->PIDs[ind] == myRank) {
808 colMapOrdering[ind].PID = -1;
810 colMapOrdering[ind].PID = ghostedCoarseNodes->PIDs[ind];
812 colMapOrdering[ind].LID = ghostedCoarseNodes->LIDs[ind];
813 colMapOrdering[ind].lexiInd = ind;
815 std::sort(colMapOrdering.begin(), colMapOrdering.end(),
817 return ((a.
PID < b.PID) || ((a.
PID == b.PID) && (a.
GID < b.GID)));
820 Array<GO> colGIDs(dofsPerNode * myGeo->lNumGhostedNodes);
821 for (LO ind = 0; ind < myGeo->lNumGhostedNodes; ++ind) {
823 ghostedCoarseNodes->colInds[colMapOrdering[ind].lexiInd] = ind;
824 for (LO dof = 0; dof < dofsPerNode; ++dof) {
825 colGIDs[dofsPerNode * ind + dof] = dofsPerNode * colMapOrdering[ind].GID + dof;
830 colGIDs.view(0, dofsPerNode *
831 myGeo->lNumCoarseNodes),
832 rowMapP->getIndexBase(),
836 colGIDs.view(0, colGIDs.size()),
837 rowMapP->getIndexBase(),
841 std::vector<size_t> strideInfo(1);
842 strideInfo[0] = dofsPerNode;
848 RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoords->getMap()->lib(),
849 myGeo->gNumCoarseNodes,
850 coarseNodesGIDs[0](),
851 fineCoords->getMap()->getIndexBase(),
852 fineCoords->getMap()->getComm());
853 RCP<const Map> coarseCoordsFineMap = MapFactory::Build(fineCoords->getMap()->lib(),
854 myGeo->gNumCoarseNodes,
855 coarseNodesGIDs[1](),
856 fineCoords->getMap()->getIndexBase(),
857 fineCoords->getMap()->getComm());
860 coarseCoordsFineMap);
862 myGeo->numDimensions);
863 coarseCoords->doImport(*fineCoords, *coarseImporter,
Xpetra::INSERT);
864 coarseCoords->replaceMap(coarseCoordsMap);
869 ghostedCoarseNodes->GIDs(),
870 fineCoords->getMap()->getIndexBase(),
871 fineCoords->getMap()->getComm());
872 RCP<const Import> ghostImporter = ImportFactory::Build(fineCoords->getMap(), ghostMap);
874 myGeo->numDimensions);
875 ghostCoords->doImport(*fineCoords, *ghostImporter,
Xpetra::INSERT);
877 P =
rcp(
new CrsMatrixWrap(rowMapP, colMapP, 0));
878 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
884 PCrs->allocateAllValues(nnzP, iaP, jaP, valP);
894 for (
int dim = 0; dim < 3; ++dim) {
895 if (dim < myGeo->numDimensions) {
896 ghostedCoords[dim] = ghostCoords->getDataNonConst(dim);
898 ghostedCoords[dim] = tmp;
907 for (
int dim = 0; dim < 3; ++dim) {
908 if (dim < myGeo->numDimensions) {
909 lFineCoords[dim] = fineCoords->getDataNonConst(dim);
911 lFineCoords[dim] = zeros->getDataNonConst(0);
916 for (
int currentIndex = 0; currentIndex < myGeo->lNumFineNodes; ++currentIndex) {
917 Array<GO> ghostedIndices(3), firstInterpolationIndices(3);
918 Array<LO> interpolationPIDs(8), interpolationLIDs(8), interpolationGIDs(8);
920 interpolationCoords[0].
resize(3);
921 GO firstInterpolationNodeIndex;
923 for (
int dim = 0; dim < 3; ++dim) {
924 interpolationCoords[0][dim] = lFineCoords[dim][currentIndex];
930 ghostedIndices[2] = currentIndex / (myGeo->lFineNodesPerDir[1] * myGeo->lFineNodesPerDir[0]);
931 LO tmp = currentIndex % (myGeo->lFineNodesPerDir[1] * myGeo->lFineNodesPerDir[0]);
932 ghostedIndices[1] = tmp / myGeo->lFineNodesPerDir[0];
933 ghostedIndices[0] = tmp % myGeo->lFineNodesPerDir[0];
934 for (
int dim = 0; dim < 3; ++dim) {
935 ghostedIndices[dim] += myGeo->offsets[dim];
944 for (
int dim = 0; dim < 3; ++dim) {
945 firstInterpolationIndices[dim] = ghostedIndices[dim] / myGeo->coarseRate[dim];
948 if (firstInterpolationIndices[dim] == myGeo->ghostedCoarseNodesPerDir[dim] - 1 && myGeo->ghostedCoarseNodesPerDir[dim] > 1) {
949 firstInterpolationIndices[dim] -= 1;
952 firstInterpolationNodeIndex =
953 firstInterpolationIndices[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + firstInterpolationIndices[1] * myGeo->ghostedCoarseNodesPerDir[0] + firstInterpolationIndices[0];
959 for (
int k = 0; k < 2; ++k) {
960 for (
int j = 0; j < 2; ++j) {
961 for (
int i = 0; i < 2; ++i) {
962 ind = k * 4 + j * 2 + i;
963 interpolationLIDs[ind] = firstInterpolationNodeIndex + k * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + j * myGeo->ghostedCoarseNodesPerDir[0] + i;
964 if (ghostedCoarseNodes->PIDs[interpolationLIDs[ind]] == rowMapP->getComm()->getRank()) {
965 interpolationPIDs[ind] = -1;
967 interpolationPIDs[ind] = ghostedCoarseNodes->PIDs[interpolationLIDs[ind]];
969 interpolationGIDs[ind] = ghostedCoarseNodes->coarseGIDs[interpolationLIDs[ind]];
971 interpolationCoords[ind + 1].
resize(3);
972 for (
int dim = 0; dim < 3; ++dim) {
973 interpolationCoords[ind + 1][dim] = ghostedCoords[dim][interpolationLIDs[ind]];
982 std::vector<double> stencil(8);
985 for (
int dim = 0; dim < 3; ++dim) {
986 firstCoarseNodeFineIndices[dim] = firstInterpolationIndices[dim] * myGeo->coarseRate[dim];
987 if ((myGeo->startIndices[dim + 3] == myGeo->gFineNodesPerDir[dim] - 1) && (ghostedIndices[dim] >=
988 (myGeo->ghostedCoarseNodesPerDir[dim] - 2) * myGeo->coarseRate[dim])) {
989 rate[dim] = myGeo->endRate[dim];
991 rate[dim] = myGeo->coarseRate[dim];
998 for (
int dim = 0; dim < 3; ++dim) {
999 if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1) {
1000 trueGhostedIndices[dim] = ghostedIndices[dim] + rate[dim];
1002 trueGhostedIndices[dim] = ghostedIndices[dim];
1005 ComputeStencil(myGeo->numDimensions, trueGhostedIndices, firstCoarseNodeFineIndices, rate,
1006 interpolationCoords, interpolationOrder, stencil);
1011 Array<LO> nzIndStencil(8), permutation(8);
1012 for (LO k = 0; k < 8; ++k) {
1015 if (interpolationOrder == 0) {
1017 for (LO k = 0; k < 8; ++k) {
1018 nzIndStencil[k] =
static_cast<LO
>(stencil[0]);
1021 stencil[nzIndStencil[0]] = 1.0;
1022 }
else if (interpolationOrder == 1) {
1023 Array<GO> currentNodeGlobalFineIndices(3);
1024 for (
int dim = 0; dim < 3; ++dim) {
1025 currentNodeGlobalFineIndices[dim] = ghostedIndices[dim] - myGeo->offsets[dim] + myGeo->startIndices[dim];
1027 if (((ghostedIndices[0] % myGeo->coarseRate[0] == 0) || currentNodeGlobalFineIndices[0] == myGeo->gFineNodesPerDir[0] - 1) && ((ghostedIndices[1] % myGeo->coarseRate[1] == 0) || currentNodeGlobalFineIndices[1] == myGeo->gFineNodesPerDir[1] - 1) && ((ghostedIndices[2] % myGeo->coarseRate[2] == 0) || currentNodeGlobalFineIndices[2] == myGeo->gFineNodesPerDir[2] - 1)) {
1028 if ((currentNodeGlobalFineIndices[0] == myGeo->gFineNodesPerDir[0] - 1) ||
1029 (ghostedIndices[0] / myGeo->coarseRate[0] == myGeo->ghostedCoarseNodesPerDir[0] - 1)) {
1030 nzIndStencil[0] += 1;
1032 if (((currentNodeGlobalFineIndices[1] == myGeo->gFineNodesPerDir[1] - 1) ||
1033 (ghostedIndices[1] / myGeo->coarseRate[1] == myGeo->ghostedCoarseNodesPerDir[1] - 1)) &&
1034 (myGeo->numDimensions > 1)) {
1035 nzIndStencil[0] += 2;
1037 if (((currentNodeGlobalFineIndices[2] == myGeo->gFineNodesPerDir[2] - 1) ||
1038 (ghostedIndices[2] / myGeo->coarseRate[2] == myGeo->ghostedCoarseNodesPerDir[2] - 1)) &&
1039 myGeo->numDimensions > 2) {
1040 nzIndStencil[0] += 4;
1043 for (LO k = 0; k < 8; ++k) {
1044 nzIndStencil[k] = nzIndStencil[0];
1048 for (LO k = 0; k < 8; ++k) {
1049 nzIndStencil[k] = k;
1061 sh_sort2(interpolationPIDs.begin(), interpolationPIDs.end(),
1062 permutation.
begin(), permutation.
end());
1065 for (LO j = 0; j < dofsPerNode; ++j) {
1066 ia[currentIndex * dofsPerNode + j + 1] = ia[currentIndex * dofsPerNode + j] + nStencil;
1067 for (LO k = 0; k < nStencil; ++k) {
1068 colInd = ghostedCoarseNodes->colInds[interpolationLIDs[nzIndStencil[permutation[k]]]];
1069 ja[ia[currentIndex * dofsPerNode + j] + k] = colInd * dofsPerNode + j;
1070 val[ia[currentIndex * dofsPerNode + j] + k] = stencil[nzIndStencil[permutation[k]]];
1073 tStencil += nStencil;
1086 PCrs->setAllValues(iaP, jaP, valP);
1087 PCrs->expertStaticFillComplete(domainMapP, rowMapP);
1565 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1567 const LO numDimensions,
const Array<GO> currentNodeIndices,
1568 const Array<GO> coarseNodeIndices,
const LO rate[3],
1570 std::vector<double>& stencil)
const {
1573 "The interpolation order can be set to 0 or 1 only.");
1575 if (interpolationOrder == 0) {
1576 ComputeConstantInterpolationStencil(numDimensions, currentNodeIndices, coarseNodeIndices,
1578 }
else if (interpolationOrder == 1) {
1579 ComputeLinearInterpolationStencil(numDimensions, coord, stencil);
1584 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1587 const Array<GO> coarseNodeIndices,
const LO rate[3],
1588 std::vector<double>& stencil)
const {
1590 if (numDimensions > 2) {
1591 if ((currentNodeIndices[2] - coarseNodeIndices[2]) > (rate[2] / 2)) {
1595 if (numDimensions > 1) {
1596 if ((currentNodeIndices[1] - coarseNodeIndices[1]) > (rate[1] / 2)) {
1600 if ((currentNodeIndices[0] - coarseNodeIndices[0]) > (rate[0] / 2)) {
1603 stencil[0] = coarseNode;
1607 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1610 std::vector<double>& stencil)
1632 LO numTerms = std::pow(2, numDimensions),
iter = 0, max_iter = 5;
1633 double functions[4][8], norm_ref = 1, norm2 = 1, tol = 1e-5;
1634 paramCoords.
size(numDimensions);
1636 while ((
iter < max_iter) && (norm2 > tol * norm_ref)) {
1639 solutionDirection.
size(numDimensions);
1640 residual.
size(numDimensions);
1644 GetInterpolationFunctions(numDimensions, paramCoords, functions);
1645 for (
LO i = 0; i < numDimensions; ++i) {
1646 residual(i) = coord[0][i];
1647 for (
LO k = 0; k < numTerms; ++k) {
1648 residual(i) -= functions[0][k] * coord[k + 1][i];
1651 norm_ref += residual(i) * residual(i);
1652 if (i == numDimensions - 1) {
1653 norm_ref = std::sqrt(norm_ref);
1657 for (
LO j = 0; j < numDimensions; ++j) {
1658 for (
LO k = 0; k < numTerms; ++k) {
1659 Jacobian(i, j) += functions[j + 1][k] * coord[k + 1][i];
1671 for (
LO i = 0; i < numDimensions; ++i) {
1672 paramCoords(i) = paramCoords(i) + solutionDirection(i);
1676 GetInterpolationFunctions(numDimensions, paramCoords, functions);
1677 for (
LO i = 0; i < numDimensions; ++i) {
1678 double tmp = coord[0][i];
1679 for (
LO k = 0; k < numTerms; ++k) {
1680 tmp -= functions[0][k] * coord[k + 1][i];
1685 norm2 = std::sqrt(norm2);
1689 for (
LO i = 0; i < 8; ++i) {
1690 stencil[i] = functions[0][i];
1695 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1699 double functions[4][8])
const {
1700 double xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
1701 if (numDimensions == 1) {
1704 }
else if (numDimensions == 2) {
1706 eta = parameters[1];
1708 }
else if (numDimensions == 3) {
1710 eta = parameters[1];
1711 zeta = parameters[2];
1715 functions[0][0] = (1.0 - xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
1716 functions[0][1] = (1.0 + xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
1717 functions[0][2] = (1.0 - xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
1718 functions[0][3] = (1.0 + xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
1719 functions[0][4] = (1.0 - xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
1720 functions[0][5] = (1.0 + xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
1721 functions[0][6] = (1.0 - xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
1722 functions[0][7] = (1.0 + xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
1724 functions[1][0] = -(1.0 - eta) * (1.0 - zeta) / denominator;
1725 functions[1][1] = (1.0 - eta) * (1.0 - zeta) / denominator;
1726 functions[1][2] = -(1.0 + eta) * (1.0 - zeta) / denominator;
1727 functions[1][3] = (1.0 + eta) * (1.0 - zeta) / denominator;
1728 functions[1][4] = -(1.0 - eta) * (1.0 + zeta) / denominator;
1729 functions[1][5] = (1.0 - eta) * (1.0 + zeta) / denominator;
1730 functions[1][6] = -(1.0 + eta) * (1.0 + zeta) / denominator;
1731 functions[1][7] = (1.0 + eta) * (1.0 + zeta) / denominator;
1733 functions[2][0] = -(1.0 - xi) * (1.0 - zeta) / denominator;
1734 functions[2][1] = -(1.0 + xi) * (1.0 - zeta) / denominator;
1735 functions[2][2] = (1.0 - xi) * (1.0 - zeta) / denominator;
1736 functions[2][3] = (1.0 + xi) * (1.0 - zeta) / denominator;
1737 functions[2][4] = -(1.0 - xi) * (1.0 + zeta) / denominator;
1738 functions[2][5] = -(1.0 + xi) * (1.0 + zeta) / denominator;
1739 functions[2][6] = (1.0 - xi) * (1.0 + zeta) / denominator;
1740 functions[2][7] = (1.0 + xi) * (1.0 + zeta) / denominator;
1742 functions[3][0] = -(1.0 - xi) * (1.0 - eta) / denominator;
1743 functions[3][1] = -(1.0 + xi) * (1.0 - eta) / denominator;
1744 functions[3][2] = -(1.0 - xi) * (1.0 + eta) / denominator;
1745 functions[3][3] = -(1.0 + xi) * (1.0 + eta) / denominator;
1746 functions[3][4] = (1.0 - xi) * (1.0 - eta) / denominator;
1747 functions[3][5] = (1.0 + xi) * (1.0 - eta) / denominator;
1748 functions[3][6] = (1.0 - xi) * (1.0 + eta) / denominator;
1749 functions[3][7] = (1.0 + xi) * (1.0 + eta) / denominator;
1753 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1759 typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1760 DT n = last1 - first1;
1765 for (DT j = 0; j < max; j++) {
1766 for (DT k = j; k >= 0; k -= m) {
1767 if (first1[first2[k + m]] >= first1[first2[k]])
1769 std::swap(first2[k + m], first2[k]);
1776 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1782 typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1783 DT n = last1 - first1;
1788 for (DT j = 0; j < max; j++) {
1789 for (DT k = j; k >= 0; k -= m) {
1790 if (first1[k + m] >= first1[k])
1792 std::swap(first1[k + m], first1[k]);
1793 std::swap(first2[k + m], first2[k]);
1800 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1804 const LO myRankIndex,
const LO pi,
const LO pj,
const LO ,
1805 const typename std::vector<std::vector<GO> >::iterator blockStart,
1806 const typename std::vector<std::vector<GO> >::iterator blockEnd,
1807 GO& myGID,
LO& myPID,
LO& myLID)
const {
1808 LO ni = -1, nj = -1, li = -1, lj = -1, lk = -1;
1809 LO myRankGuess = myRankIndex;
1811 if (i == 0 && myGeo->ghostInterface[0]) {
1813 }
else if ((i == myGeo->ghostedCoarseNodesPerDir[0] - 1) && myGeo->ghostInterface[1]) {
1816 if (j == 0 && myGeo->ghostInterface[2]) {
1818 }
else if ((j == myGeo->ghostedCoarseNodesPerDir[1] - 1) && myGeo->ghostInterface[3]) {
1821 if (k == 0 && myGeo->ghostInterface[4]) {
1822 myRankGuess -= pj * pi;
1823 }
else if ((k == myGeo->ghostedCoarseNodesPerDir[2] - 1) && myGeo->ghostInterface[5]) {
1824 myRankGuess += pj * pi;
1826 if (coarseNodeFineIndices[0] >= myGeo->meshData[myRankGuess][3] && coarseNodeFineIndices[0] <= myGeo->meshData[myRankGuess][4] && coarseNodeFineIndices[1] >= myGeo->meshData[myRankGuess][5] && coarseNodeFineIndices[1] <= myGeo->meshData[myRankGuess][6] && coarseNodeFineIndices[2] >= myGeo->meshData[myRankGuess][7] && coarseNodeFineIndices[2] <= myGeo->meshData[myRankGuess][8]) {
1827 myPID = myGeo->meshData[myRankGuess][0];
1828 ni = myGeo->meshData[myRankGuess][4] - myGeo->meshData[myRankGuess][3] + 1;
1829 nj = myGeo->meshData[myRankGuess][6] - myGeo->meshData[myRankGuess][5] + 1;
1830 li = coarseNodeFineIndices[0] - myGeo->meshData[myRankGuess][3];
1831 lj = coarseNodeFineIndices[1] - myGeo->meshData[myRankGuess][5];
1832 lk = coarseNodeFineIndices[2] - myGeo->meshData[myRankGuess][7];
1833 myLID = lk * nj * ni + lj * ni + li;
1834 myGID = myGeo->meshData[myRankGuess][9] + myLID;
1838 auto nodeRank = std::find_if(blockStart, blockEnd,
1839 [coarseNodeFineIndices](
const std::vector<GO>& vec) {
1840 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]) {
1846 myPID = (*nodeRank)[0];
1847 ni = (*nodeRank)[4] - (*nodeRank)[3] + 1;
1848 nj = (*nodeRank)[6] - (*nodeRank)[5] + 1;
1849 li = coarseNodeFineIndices[0] - (*nodeRank)[3];
1850 lj = coarseNodeFineIndices[1] - (*nodeRank)[5];
1851 lk = coarseNodeFineIndices[2] - (*nodeRank)[7];
1852 myLID = lk * nj * ni + lj * ni + li;
1853 myGID = (*nodeRank)[9] + myLID;
1859 #define MUELU_GENERALGEOMETRICPFACTORY_SHORT
1860 #endif // MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
void ComputeLinearInterpolationStencil(const LO numDimension, const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, std::vector< double > &stencil) const
void GetCoarsePoints(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
void MakeGeneralGeometricP(RCP< GeometricData > myGeo, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &fCoords, const LO nnzP, const LO dofsPerNode, RCP< const Map > &stridedDomainMapP, RCP< Matrix > &Amat, RCP< Matrix > &P, RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &cCoords, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > coarseNodesGIDs, int interpolationOrder) const
T & get(const std::string &name, T def_value)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static const NoFactory * get()
void resize(const size_type n, const T &val=T())
void MeshLayoutInterface(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
void factorWithEquilibration(bool flag)
void sh_sort2(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
void ComputeStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, const int interpolationOrder, std::vector< double > &stencil) const
void GetGIDLocalLexicographic(const GO i, const GO j, const GO k, const Array< LO > coarseNodeFineIndices, const RCP< GeometricData > myGeo, const LO myRankIndex, const LO pi, const LO pj, const LO pk, const typename std::vector< std::vector< GO > >::iterator blockStart, const typename std::vector< std::vector< GO > >::iterator blockEnd, GO &myGID, LO &myPID, LO &myLID) const
void resize(size_type new_size, const value_type &x=value_type())
void GetInterpolationFunctions(const LO numDimension, const Teuchos::SerialDenseVector< LO, double > parameters, double functions[4][8]) const
void ComputeConstantInterpolationStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], std::vector< double > &stencil) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void sh_sort_permute(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
int size(OrdinalType length_in)
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
int GetLevelID() const
Return level number.
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
std::vector< T >::iterator iterator
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)