140 Level& coarseLevel)
const{
144 const ParameterList& pL = GetParameterList();
147 RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel,
"A");
148 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
149 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > coordinates =
150 Get< RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > >(fineLevel,
"Coordinates");
151 LO numDimensions = coordinates->getNumVectors();
152 LO BlkSize = A->GetFixedBlockSize();
155 Array<GO> gFineNodesPerDir(3);
156 Array<LO> lFineNodesPerDir(3);
163 gFineNodesPerDir = Get<Array<GO> >(fineLevel,
"gNodesPerDim");
166 lFineNodesPerDir = Get<Array<LO> >(fineLevel,
"lNodesPerDim");
168 for(LO i = 0; i < 3; ++i) {
169 if(gFineNodesPerDir[i] == 0) {
170 GetOStream(
Runtime0) <<
"gNodesPerDim in direction " << i <<
" is set to 1 from 0"
172 gFineNodesPerDir[i] = 1;
174 if(lFineNodesPerDir[i] == 0) {
175 GetOStream(
Runtime0) <<
"lNodesPerDim in direction " << i <<
" is set to 1 from 0"
177 lFineNodesPerDir[i] = 1;
180 GO gNumFineNodes = gFineNodesPerDir[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
181 LO lNumFineNodes = lFineNodesPerDir[2]*lFineNodesPerDir[1]*lFineNodesPerDir[0];
184 std::string coarsenRate = pL.get<std::string>(
"Coarsen");
185 Array<LO> coarseRate(3);
187 Teuchos::Array<LO> crates;
189 crates = Teuchos::fromStringToArray<LO>(coarsenRate);
190 }
catch(
const Teuchos::InvalidArrayStringRepresentation& e) {
191 GetOStream(
Errors,-1) <<
" *** Coarsen must be a string convertible into an array! *** "
195 TEUCHOS_TEST_FOR_EXCEPTION((crates.size() > 1) && (crates.size() < numDimensions),
197 "Coarsen must have at least as many components as the number of"
198 " spatial dimensions in the problem.");
199 for(LO i = 0; i < 3; ++i) {
200 if(i < numDimensions) {
201 if(crates.size() == 1) {
202 coarseRate[i] = crates[0];
203 }
else if(i < crates.size()) {
204 coarseRate[i] = crates[i];
206 GetOStream(
Errors,-1) <<
" *** Coarsen must be at least as long as the number of"
207 " spatial dimensions! *** " << std::endl;
209 " spatial dimensions! *** \n");
218 const std::string stencilType = pL.get<std::string>(
"stencil type");
219 if(stencilType !=
"full" && stencilType !=
"reduced") {
220 GetOStream(
Errors,-1) <<
" *** stencil type must be set to: full or reduced, any other value "
221 "is trated as an error! *** " << std::endl;
226 const std::string blockStrategy = pL.get<std::string>(
"block strategy");
227 if(blockStrategy !=
"coupled" && blockStrategy !=
"uncoupled") {
228 GetOStream(
Errors,-1) <<
" *** block strategy must be set to: coupled or uncoupled, any other "
229 "value is trated as an error! *** " << std::endl;
233 GO gNumCoarseNodes = 0;
234 LO lNumCoarseNodes = 0;
235 Array<GO> gIndices(3), gCoarseNodesPerDir(3), ghostGIDs, coarseNodesGIDs, colGIDs;
236 Array<LO> myOffset(3), lCoarseNodesPerDir(3), glCoarseNodesPerDir(3), endRate(3);
237 Array<bool> ghostInterface(6);
238 Array<int> boundaryFlags(3);
239 ArrayRCP<Array<typename Teuchos::ScalarTraits<Scalar>::magnitudeType> > coarseNodes(numDimensions);
240 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType> > fineNodes(numDimensions);
241 for(LO dim = 0; dim < numDimensions; ++dim) {fineNodes[dim] = coordinates->getData(dim)();}
244 RCP<NodesIDs> ghostedCoarseNodes = rcp(
new NodesIDs{});
246 GetGeometricData(coordinates, coarseRate, gFineNodesPerDir, lFineNodesPerDir, BlkSize,
247 gIndices, myOffset, ghostInterface, endRate, gCoarseNodesPerDir,
248 lCoarseNodesPerDir, glCoarseNodesPerDir, ghostGIDs, coarseNodesGIDs, colGIDs,
249 gNumCoarseNodes, lNumCoarseNodes, coarseNodes, boundaryFlags,
253 Xpetra::UnderlyingLib lib = coordinates->getMap()->lib();
254 RCP<const Map> coarseCoordsMap = MapFactory::Build (lib,
256 coarseNodesGIDs.view(0, lNumCoarseNodes),
257 coordinates->getMap()->getIndexBase(),
258 coordinates->getMap()->getComm());
259 Array<ArrayView<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType> > coarseCoords(numDimensions);
260 for(LO dim = 0; dim < numDimensions; ++dim) {
261 coarseCoords[dim] = coarseNodes[dim]();
263 RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> > coarseCoordinates =
264 Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO>::Build(coarseCoordsMap, coarseCoords(),
297 Array<GO> ghostRowGIDs, ghostColGIDs, nodeSteps(3);
299 nodeSteps[1] = gFineNodesPerDir[0];
300 nodeSteps[2] = gFineNodesPerDir[0]*gFineNodesPerDir[1];
301 Array<LO> glFineNodesPerDir(3);
302 GO startingGID = A->getRowMap()->getMinGlobalIndex();
303 for(LO dim = 0; dim < 3; ++dim) {
304 LO numCoarseNodes = 0;
305 if(dim < numDimensions) {
306 startingGID -= myOffset[dim]*nodeSteps[dim];
307 numCoarseNodes = lCoarseNodesPerDir[dim];
308 if(ghostInterface[2*dim]) {++numCoarseNodes;}
309 if(ghostInterface[2*dim+1]) {++numCoarseNodes;}
310 if(gIndices[dim] + lFineNodesPerDir[dim] == gFineNodesPerDir[dim]) {
311 glFineNodesPerDir[dim] = (numCoarseNodes - 2)*coarseRate[dim] + endRate[dim] + 1;
313 glFineNodesPerDir[dim] = (numCoarseNodes - 1)*coarseRate[dim] + 1;
316 glFineNodesPerDir[dim] = 1;
319 ghostRowGIDs.resize(glFineNodesPerDir[0]*glFineNodesPerDir[1]*glFineNodesPerDir[2]*BlkSize);
320 for(LO k = 0; k < glFineNodesPerDir[2]; ++k) {
321 for(LO j = 0; j < glFineNodesPerDir[1]; ++j) {
322 for(LO i = 0; i < glFineNodesPerDir[0]; ++i) {
323 for(LO l = 0; l < BlkSize; ++l) {
324 ghostRowGIDs[(k*glFineNodesPerDir[1]*glFineNodesPerDir[0]
325 + j*glFineNodesPerDir[0] + i)*BlkSize + l] = startingGID
326 + (k*gFineNodesPerDir[1]*gFineNodesPerDir[0] + j*gFineNodesPerDir[0] + i)*BlkSize + l;
333 Array<GO> startingGlobalIndices(numDimensions), dimStride(numDimensions);
334 Array<GO> startingColIndices(numDimensions), finishingColIndices(numDimensions);
336 Array<LO> colRange(numDimensions);
338 for(
int dim = 1; dim < numDimensions; ++dim) {
339 dimStride[dim] = dimStride[dim - 1]*gFineNodesPerDir[dim - 1];
342 GO tmp = startingGID;
343 for(
int dim = numDimensions; dim > 0; --dim) {
344 startingGlobalIndices[dim - 1] = tmp / dimStride[dim - 1];
345 tmp = tmp % dimStride[dim - 1];
347 if(startingGlobalIndices[dim - 1] > 0) {
348 startingColIndices[dim - 1] = startingGlobalIndices[dim - 1] - 1;
350 if(startingGlobalIndices[dim - 1] + glFineNodesPerDir[dim - 1] < gFineNodesPerDir[dim - 1]){
351 finishingColIndices[dim - 1] = startingGlobalIndices[dim - 1]
352 + glFineNodesPerDir[dim - 1];
354 finishingColIndices[dim - 1] = startingGlobalIndices[dim - 1]
355 + glFineNodesPerDir[dim - 1] - 1;
357 colRange[dim - 1] = finishingColIndices[dim - 1] - startingColIndices[dim - 1] + 1;
358 colMinGID += startingColIndices[dim - 1]*dimStride[dim - 1];
361 ghostColGIDs.resize(colRange[0]*colRange[1]*colRange[2]*BlkSize);
362 for(LO k = 0; k < colRange[2]; ++k) {
363 for(LO j = 0; j < colRange[1]; ++j) {
364 for(LO i = 0; i < colRange[0]; ++i) {
365 for(LO l = 0; l < BlkSize; ++l) {
366 ghostColGIDs[(k*colRange[1]*colRange[0] + j*colRange[0] + i)*BlkSize + l] = colMinGID
367 + (k*gFineNodesPerDir[1]*gFineNodesPerDir[0] + j*gFineNodesPerDir[0] + i)*BlkSize + l;
373 RCP<const Map> ghostedRowMap = Xpetra::MapFactory<LO,GO,NO>::Build(A->getRowMap()->lib(),
374 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
376 A->getRowMap()->getIndexBase(),
377 A->getRowMap()->getComm());
378 RCP<const Map> ghostedColMap = Xpetra::MapFactory<LO,GO,NO>::Build(A->getRowMap()->lib(),
379 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
381 A->getRowMap()->getIndexBase(),
382 A->getRowMap()->getComm());
383 RCP<const Import> ghostImporter = Xpetra::ImportFactory<LO,GO,NO>::Build(A->getRowMap(),
385 RCP<const Matrix> Aghost = Xpetra::MatrixFactory<SC,LO,GO,NO>::Build(A, *ghostImporter,
390 RCP<const Map> rowMapP = A->getDomainMap();
392 RCP<const Map> domainMapP;
394 RCP<const Map> colMapP;
405 LO lNumGhostedNodes = ghostedCoarseNodes->GIDs.size();
407 std::vector<NodeID> colMapOrdering(lNumGhostedNodes);
408 for(LO ind = 0; ind < lNumGhostedNodes; ++ind) {
409 colMapOrdering[ind].GID = ghostedCoarseNodes->GIDs[ind];
410 if(ghostedCoarseNodes->PIDs[ind] == rowMapP->getComm()->getRank()) {
411 colMapOrdering[ind].PID = -1;
413 colMapOrdering[ind].PID = ghostedCoarseNodes->PIDs[ind];
415 colMapOrdering[ind].LID = ghostedCoarseNodes->LIDs[ind];
416 colMapOrdering[ind].lexiInd = ind;
418 std::sort(colMapOrdering.begin(), colMapOrdering.end(),
420 return ( (a.PID < b.PID) || ((a.PID == b.PID) && (a.GID < b.GID)) );
423 colGIDs.resize(BlkSize*lNumGhostedNodes);
424 for(LO ind = 0; ind < lNumGhostedNodes; ++ind) {
426 ghostedCoarseNodes->colInds[colMapOrdering[ind].lexiInd] = ind;
427 for(LO dof = 0; dof < BlkSize; ++dof) {
428 colGIDs[BlkSize*ind + dof] = BlkSize*colMapOrdering[ind].GID + dof;
431 domainMapP = Xpetra::MapFactory<LO,GO,NO>::Build(rowMapP->lib(),
432 BlkSize*gNumCoarseNodes,
433 colGIDs.view(0,BlkSize*lNumCoarseNodes),
434 rowMapP->getIndexBase(),
436 colMapP = Xpetra::MapFactory<LO,GO,NO>::Build(rowMapP->lib(),
437 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
438 colGIDs.view(0, colGIDs.size()),
439 rowMapP->getIndexBase(),
443 std::vector<size_t> strideInfo(1);
444 strideInfo[0] = BlkSize;
445 RCP<const Map> stridedDomainMapP = Xpetra::StridedMapFactory<LO,GO,NO>::Build(domainMapP,
451 gnnzP += gNumCoarseNodes;
452 lnnzP += lNumCoarseNodes;
454 gnnzP += (gNumFineNodes - gNumCoarseNodes)*std::pow(2, numDimensions);
455 lnnzP += (lNumFineNodes - lNumCoarseNodes)*std::pow(2, numDimensions);
457 gnnzP = gnnzP*BlkSize;
458 lnnzP = lnnzP*BlkSize;
462 P = rcp(
new CrsMatrixWrap(rowMapP, colMapP, 0));
463 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
465 ArrayRCP<size_t> iaP;
469 PCrs->allocateAllValues(lnnzP, iaP, jaP, valP);
471 ArrayView<size_t> ia = iaP();
472 ArrayView<LO> ja = jaP();
473 ArrayView<SC> val = valP();
477 LO numCoarseElements = 1;
478 Array<LO> lCoarseElementsPerDir(3);
479 for(LO dim = 0; dim < numDimensions; ++dim) {
480 lCoarseElementsPerDir[dim] = lCoarseNodesPerDir[dim];
481 if(ghostInterface[2*dim]) {++lCoarseElementsPerDir[dim];}
482 if(!ghostInterface[2*dim+1]) {--lCoarseElementsPerDir[dim];}
483 numCoarseElements = numCoarseElements*lCoarseElementsPerDir[dim];
486 for(LO dim = numDimensions; dim < 3; ++dim) {
487 lCoarseElementsPerDir[dim] = 1;
491 Array<int> elementFlags(3);
492 Array<LO> elemInds(3), elementNodesPerDir(3), glElementRefTuple(3);
493 Array<LO> glElementRefTupleCG(3), glElementCoarseNodeCG(8);
494 const int numCoarseNodesInElement = std::pow(2, numDimensions);
495 const int nnzPerCoarseNode = (blockStrategy ==
"coupled") ? BlkSize : 1;
496 const int numRowsPerPoint = BlkSize;
497 for(elemInds[2] = 0; elemInds[2] < lCoarseElementsPerDir[2]; ++elemInds[2]) {
498 for(elemInds[1] = 0; elemInds[1] < lCoarseElementsPerDir[1]; ++elemInds[1]) {
499 for(elemInds[0] = 0; elemInds[0] < lCoarseElementsPerDir[0]; ++elemInds[0]) {
500 elementFlags[0] = 0; elementFlags[1] = 0; elementFlags[2] = 0;
501 for(
int dim = 0; dim < 3; ++dim) {
503 if(elemInds[dim] == 0 && elemInds[dim] == lCoarseElementsPerDir[dim] - 1) {
504 elementFlags[dim] = boundaryFlags[dim];
505 }
else if(elemInds[dim] == 0 && (boundaryFlags[dim] == 1 || boundaryFlags[dim] == 3)) {
506 elementFlags[dim] += 1;
507 }
else if((elemInds[dim] == lCoarseElementsPerDir[dim] - 1)
508 && (boundaryFlags[dim] == 2 || boundaryFlags[dim] == 3)) {
509 elementFlags[dim] += 2;
511 elementFlags[dim] = 0;
515 if(dim < numDimensions) {
516 if((elemInds[dim] == lCoarseElementsPerDir[dim])
517 && (gIndices[dim] + lFineNodesPerDir[dim] == gFineNodesPerDir[dim])) {
518 elementNodesPerDir[dim] = endRate[dim] + 1;
520 elementNodesPerDir[dim] = coarseRate[dim] + 1;
523 elementNodesPerDir[dim] = 1;
527 glElementRefTuple[dim] = elemInds[dim]*coarseRate[dim];
528 glElementRefTupleCG[dim] = elemInds[dim];
533 for(
typename Array<LO>::size_type elem = 0; elem < glElementCoarseNodeCG.size(); ++elem) {
534 glElementCoarseNodeCG[elem]
535 = glElementRefTupleCG[2]*glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0]
536 + glElementRefTupleCG[1]*glCoarseNodesPerDir[0] + glElementRefTupleCG[0];
538 glElementCoarseNodeCG[4] += glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0];
539 glElementCoarseNodeCG[5] += glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0];
540 glElementCoarseNodeCG[6] += glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0];
541 glElementCoarseNodeCG[7] += glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0];
543 glElementCoarseNodeCG[2] += glCoarseNodesPerDir[0];
544 glElementCoarseNodeCG[3] += glCoarseNodesPerDir[0];
545 glElementCoarseNodeCG[6] += glCoarseNodesPerDir[0];
546 glElementCoarseNodeCG[7] += glCoarseNodesPerDir[0];
548 glElementCoarseNodeCG[1] += 1;
549 glElementCoarseNodeCG[3] += 1;
550 glElementCoarseNodeCG[5] += 1;
551 glElementCoarseNodeCG[7] += 1;
553 LO numNodesInElement = elementNodesPerDir[0]*elementNodesPerDir[1]*elementNodesPerDir[2];
558 Teuchos::SerialDenseMatrix<LO,SC> Pi, Pf, Pe;
559 Array<LO> dofType(numNodesInElement*BlkSize), lDofInd(numNodesInElement*BlkSize);
560 ComputeLocalEntries(Aghost, coarseRate, endRate, BlkSize, elemInds, lCoarseElementsPerDir,
561 numDimensions, glFineNodesPerDir, gFineNodesPerDir, gIndices,
562 lCoarseNodesPerDir, ghostInterface, elementFlags, stencilType,
563 blockStrategy, elementNodesPerDir, numNodesInElement, colGIDs,
564 Pi, Pf, Pe, dofType, lDofInd);
568 Array<LO> lNodeLIDs(numNodesInElement);
570 Array<LO> lNodeTuple(3), nodeInd(3);
571 for(nodeInd[2] = 0; nodeInd[2] < elementNodesPerDir[2]; ++nodeInd[2]) {
572 for(nodeInd[1] = 0; nodeInd[1] < elementNodesPerDir[1]; ++nodeInd[1]) {
573 for(nodeInd[0] = 0; nodeInd[0] < elementNodesPerDir[0]; ++nodeInd[0]) {
574 int stencilLength = 0;
575 if((nodeInd[0] == 0 || nodeInd[0] == elementNodesPerDir[0] - 1) &&
576 (nodeInd[1] == 0 || nodeInd[1] == elementNodesPerDir[1] - 1) &&
577 (nodeInd[2] == 0 || nodeInd[2] == elementNodesPerDir[2] - 1)) {
580 stencilLength = std::pow(2, numDimensions);
582 LO nodeElementInd = nodeInd[2]*elementNodesPerDir[1]*elementNodesPerDir[1]
583 + nodeInd[1]*elementNodesPerDir[0] + nodeInd[0];
584 for(
int dim = 0; dim < 3; ++dim) {
585 lNodeTuple[dim] = glElementRefTuple[dim] - myOffset[dim] + nodeInd[dim];
587 if(lNodeTuple[0] < 0 || lNodeTuple[1] < 0 || lNodeTuple[2] < 0 ||
588 lNodeTuple[0] > lFineNodesPerDir[0] - 1 ||
589 lNodeTuple[1] > lFineNodesPerDir[1] - 1 ||
590 lNodeTuple[2] > lFineNodesPerDir[2] - 1) {
593 lNodeLIDs[nodeElementInd] = -1;
594 }
else if ((nodeInd[0] == 0 && elemInds[0] !=0) ||
595 (nodeInd[1] == 0 && elemInds[1] !=0) ||
596 (nodeInd[2] == 0 && elemInds[2] !=0)) {
600 lNodeLIDs[nodeElementInd] = -2;
606 lNodeLIDs[nodeElementInd] = lNodeTuple[2]*lFineNodesPerDir[1]*lFineNodesPerDir[0]
607 + lNodeTuple[1]*lFineNodesPerDir[0] + lNodeTuple[0];
613 Array<LO> refCoarsePointTuple(3);
614 for(
int dim = 2; dim > -1; --dim) {
616 refCoarsePointTuple[dim] = (lNodeTuple[dim] + myOffset[dim]) / coarseRate[dim];
617 if(myOffset[dim] == 0) {
618 ++refCoarsePointTuple[dim];
622 refCoarsePointTuple[dim] =
623 std::ceil(
static_cast<double>(lNodeTuple[dim] + myOffset[dim])
626 if((lNodeTuple[dim] + myOffset[dim]) % coarseRate[dim] > 0) {
break;}
628 const LO numCoarsePoints = refCoarsePointTuple[0]
629 + refCoarsePointTuple[1]*lCoarseNodesPerDir[0]
630 + refCoarsePointTuple[2]*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0];
631 const LO numFinePoints = lNodeLIDs[nodeElementInd] + 1;
635 size_t rowPtr = (numFinePoints - numCoarsePoints)*numRowsPerPoint
636 *numCoarseNodesInElement*nnzPerCoarseNode + numCoarsePoints*numRowsPerPoint;
637 if(dofType[nodeElementInd*BlkSize] == 0) {
639 rowPtr = rowPtr - numRowsPerPoint;
641 rowPtr = rowPtr - numRowsPerPoint*numCoarseNodesInElement*nnzPerCoarseNode;
644 for(
int dof = 0; dof < BlkSize; ++dof) {
648 switch(dofType[nodeElementInd*BlkSize + dof]) {
650 if(nodeInd[2] == elementNodesPerDir[2] - 1) {cornerInd += 4;}
651 if(nodeInd[1] == elementNodesPerDir[1] - 1) {cornerInd += 2;}
652 if(nodeInd[0] == elementNodesPerDir[0] - 1) {cornerInd += 1;}
653 ia[lNodeLIDs[nodeElementInd]*BlkSize + dof + 1] = rowPtr + dof + 1;
654 ja[rowPtr + dof] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[cornerInd]]*BlkSize + dof;
655 val[rowPtr + dof] = 1.0;
659 ia[lNodeLIDs[nodeElementInd]*BlkSize + dof + 1]
660 = rowPtr + (dof + 1)*numCoarseNodesInElement*nnzPerCoarseNode;
661 for(
int ind1 = 0; ind1 < stencilLength; ++ind1) {
662 if(blockStrategy ==
"coupled") {
663 for(
int ind2 = 0; ind2 < BlkSize; ++ind2) {
664 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
665 + ind1*BlkSize + ind2;
666 ja[lRowPtr] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + ind2;
667 val[lRowPtr] = Pe(lDofInd[nodeElementInd*BlkSize + dof],
668 ind1*BlkSize + ind2);
670 }
else if(blockStrategy ==
"uncoupled") {
671 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
673 ja[rowPtr + dof] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + dof;
674 val[lRowPtr] = Pe(lDofInd[nodeElementInd*BlkSize + dof],
681 ia[lNodeLIDs[nodeElementInd]*BlkSize + dof + 1]
682 = rowPtr + (dof + 1)*numCoarseNodesInElement*nnzPerCoarseNode;
683 for(
int ind1 = 0; ind1 < stencilLength; ++ind1) {
684 if(blockStrategy ==
"coupled") {
685 for(
int ind2 = 0; ind2 < BlkSize; ++ind2) {
686 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
687 + ind1*BlkSize + ind2;
688 ja [lRowPtr] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + ind2;
689 val[lRowPtr] = Pf(lDofInd[nodeElementInd*BlkSize + dof],
690 ind1*BlkSize + ind2);
692 }
else if(blockStrategy ==
"uncoupled") {
693 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
696 ja [lRowPtr] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + dof;
697 val[lRowPtr] = Pf(lDofInd[nodeElementInd*BlkSize + dof],
704 ia[lNodeLIDs[nodeElementInd]*BlkSize + dof + 1]
705 = rowPtr + (dof + 1)*numCoarseNodesInElement*nnzPerCoarseNode;
706 for(
int ind1 = 0; ind1 < stencilLength; ++ind1) {
707 if(blockStrategy ==
"coupled") {
708 for(
int ind2 = 0; ind2 < BlkSize; ++ind2) {
709 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
710 + ind1*BlkSize + ind2;
711 ja [lRowPtr] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + ind2;
712 val[lRowPtr] = Pi(lDofInd[nodeElementInd*BlkSize + dof],
713 ind1*BlkSize + ind2);
715 }
else if(blockStrategy ==
"uncoupled") {
716 size_t lRowPtr = rowPtr + dof*numCoarseNodesInElement*nnzPerCoarseNode
718 ja [lRowPtr] = ghostedCoarseNodes->colInds[glElementCoarseNodeCG[ind1]]*BlkSize + dof;
719 val[lRowPtr] = Pi(lDofInd[nodeElementInd*BlkSize + dof],
736 Xpetra::CrsMatrixUtils<SC,LO,GO,NO>::sortCrsEntries(ia, ja, val, rowMapP->lib());
739 PCrs->setAllValues(iaP, jaP, valP);
740 PCrs->expertStaticFillComplete(domainMapP,rowMapP);
743 if (A->IsView(
"stridedMaps") ==
true) {
744 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMapP);
746 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMapP);
750 Set(coarseLevel,
"P", P);
751 Set(coarseLevel,
"coarseCoordinates", coarseCoordinates);
752 Set<Array<GO> >(coarseLevel,
"gCoarseNodesPerDim", gCoarseNodesPerDir);
753 Set<Array<LO> >(coarseLevel,
"lCoarseNodesPerDim", lCoarseNodesPerDir);
759 GetGeometricData(RCP<Xpetra::MultiVector<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType,LO,GO,NO> >& coordinates,
760 const Array<LO> coarseRate,
const Array<GO> gFineNodesPerDir,
761 const Array<LO> lFineNodesPerDir,
const LO BlkSize, Array<GO>& gIndices,
762 Array<LO>& myOffset, Array<bool>& ghostInterface, Array<LO>& endRate,
763 Array<GO>& gCoarseNodesPerDir, Array<LO>& lCoarseNodesPerDir,
764 Array<LO>& glCoarseNodesPerDir, Array<GO>& ghostGIDs, Array<GO>& coarseNodesGIDs,
765 Array<GO>& colGIDs, GO& gNumCoarseNodes, LO& lNumCoarseNodes,
766 ArrayRCP<Array<
typename Teuchos::ScalarTraits<Scalar>::magnitudeType> > coarseNodes, Array<int>& boundaryFlags,
767 RCP<NodesIDs> ghostedCoarseNodes)
const {
774 RCP<const Map> coordinatesMap = coordinates->getMap();
775 LO numDimensions = coordinates->getNumVectors();
786 GO minGlobalIndex = coordinatesMap->getMinGlobalIndex();
789 gIndices[2] = minGlobalIndex / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
790 tmp = minGlobalIndex % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
791 gIndices[1] = tmp / gFineNodesPerDir[0];
792 gIndices[0] = tmp % gFineNodesPerDir[0];
794 myOffset[2] = gIndices[2] % coarseRate[2];
795 myOffset[1] = gIndices[1] % coarseRate[1];
796 myOffset[0] = gIndices[0] % coarseRate[0];
799 for(
int dim = 0; dim < 3; ++dim) {
800 if(gIndices[dim] == 0) {
801 boundaryFlags[dim] += 1;
803 if(gIndices[dim] + lFineNodesPerDir[dim] == gFineNodesPerDir[dim]) {
804 boundaryFlags[dim] += 2;
809 for(LO i=0; i < numDimensions; ++i) {
810 if((gIndices[i] != 0) && (gIndices[i] % coarseRate[i] > 0)) {
811 ghostInterface[2*i] =
true;
813 if(((gIndices[i] + lFineNodesPerDir[i]) != gFineNodesPerDir[i])
814 && ((gIndices[i] + lFineNodesPerDir[i] - 1) % coarseRate[i] > 0)) {
815 ghostInterface[2*i + 1] =
true;
819 for(LO i = 0; i < 3; ++i) {
820 if(i < numDimensions) {
821 lCoarseNodesPerDir[i] = (lFineNodesPerDir[i] + myOffset[i] - 1) / coarseRate[i];
822 if(myOffset[i] == 0) { ++lCoarseNodesPerDir[i]; }
823 gCoarseNodesPerDir[i] = (gFineNodesPerDir[i] - 1) / coarseRate[i];
824 endRate[i] = (gFineNodesPerDir[i] - 1) % coarseRate[i];
825 if(endRate[i] == 0) {
826 ++gCoarseNodesPerDir[i];
827 endRate[i] = coarseRate[i];
833 gCoarseNodesPerDir[i] = 1;
834 lCoarseNodesPerDir[i] = 1;
839 gNumCoarseNodes = gCoarseNodesPerDir[0]*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[2];
840 lNumCoarseNodes = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];
843 LO lNumGhostNodes = 0;
844 Array<GO> startGhostedCoarseNode(3);
846 const int complementaryIndices[3][2] = {{1,2}, {0,2}, {0,1}};
848 for(LO i = 0; i < 3; ++i) {
852 if((gIndices[i] == gFineNodesPerDir[i] - 1) && (gIndices[i] % coarseRate[i] == 0)) {
853 startGhostedCoarseNode[i] = gIndices[i] / coarseRate[i] - 1;
855 startGhostedCoarseNode[i] = gIndices[i] / coarseRate[i];
858 glCoarseNodesPerDir[i] = lCoarseNodesPerDir[i];
861 if(ghostInterface[2*i] || ghostInterface[2*i+1]) {
862 ++glCoarseNodesPerDir[i];
863 if(i == 0) {tmp = lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];}
864 if(i == 1) {tmp = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[2];}
865 if(i == 2) {tmp = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1];}
868 if(ghostInterface[2*i] && ghostInterface[2*i+1]) {
869 ++glCoarseNodesPerDir[i];
872 lNumGhostNodes += tmp;
875 for(LO j = 0; j < 2; ++j) {
876 for(LO k = 0; k < 2; ++k) {
878 if(ghostInterface[2*complementaryIndices[i][0]+j]
879 && ghostInterface[2*complementaryIndices[i][1]+k]) {
880 lNumGhostNodes += lCoarseNodesPerDir[i];
883 if(ghostInterface[2*i] && (i == 0)) { lNumGhostNodes += 1; }
884 if(ghostInterface[2*i+1] && (i == 0)) { lNumGhostNodes += 1; }
892 const LO lNumGhostedNodes = glCoarseNodesPerDir[0]*glCoarseNodesPerDir[1]
893 *glCoarseNodesPerDir[2];
894 ghostedCoarseNodes->PIDs.resize(lNumGhostedNodes);
895 ghostedCoarseNodes->LIDs.resize(lNumGhostedNodes);
896 ghostedCoarseNodes->GIDs.resize(lNumGhostedNodes);
897 ghostedCoarseNodes->coarseGIDs.resize(lNumGhostedNodes);
898 ghostedCoarseNodes->colInds.resize(lNumGhostedNodes);
901 Array<LO> coarseNodeCoarseIndices(3), coarseNodeFineIndices(3), ijk(3);
902 LO currentIndex = -1;
903 for(ijk[2] = 0; ijk[2] < glCoarseNodesPerDir[2]; ++ijk[2]) {
904 for(ijk[1] = 0; ijk[1] < glCoarseNodesPerDir[1]; ++ijk[1]) {
905 for(ijk[0] = 0; ijk[0] < glCoarseNodesPerDir[0]; ++ijk[0]) {
906 currentIndex = ijk[2]*glCoarseNodesPerDir[1]*glCoarseNodesPerDir[0]
907 + ijk[1]*glCoarseNodesPerDir[0] + ijk[0];
908 coarseNodeCoarseIndices[0] = startGhostedCoarseNode[0] + ijk[0];
909 coarseNodeFineIndices[0] = coarseNodeCoarseIndices[0]*coarseRate[0];
910 if(coarseNodeFineIndices[0] > gFineNodesPerDir[0] - 1) {
911 coarseNodeFineIndices[0] = gFineNodesPerDir[0] - 1;
913 coarseNodeCoarseIndices[1] = startGhostedCoarseNode[1] + ijk[1];
914 coarseNodeFineIndices[1] = coarseNodeCoarseIndices[1]*coarseRate[1];
915 if(coarseNodeFineIndices[1] > gFineNodesPerDir[1] - 1) {
916 coarseNodeFineIndices[1] = gFineNodesPerDir[1] - 1;
918 coarseNodeCoarseIndices[2] = startGhostedCoarseNode[2] + ijk[2];
919 coarseNodeFineIndices[2] = coarseNodeCoarseIndices[2]*coarseRate[2];
920 if(coarseNodeFineIndices[2] > gFineNodesPerDir[2] - 1) {
921 coarseNodeFineIndices[2] = gFineNodesPerDir[2] - 1;
923 GO myGID = 0, myCoarseGID = -1;
925 factor[2] = gFineNodesPerDir[1]*gFineNodesPerDir[0];
926 factor[1] = gFineNodesPerDir[0];
928 for(
int dim = 0; dim < 3; ++dim) {
929 if(dim < numDimensions) {
930 if(gIndices[dim] - myOffset[dim] + ijk[dim]*coarseRate[dim]
931 < gFineNodesPerDir[dim] - 1) {
932 myGID += (gIndices[dim] - myOffset[dim]
933 + ijk[dim]*coarseRate[dim])*factor[dim];
935 myGID += (gIndices[dim] - myOffset[dim]
936 + (ijk[dim] - 1)*coarseRate[dim] + endRate[dim])*factor[dim];
940 myCoarseGID = coarseNodeCoarseIndices[0]
941 + coarseNodeCoarseIndices[1]*gCoarseNodesPerDir[0]
942 + coarseNodeCoarseIndices[2]*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0];
943 ghostedCoarseNodes->GIDs[currentIndex] = myGID;
944 ghostedCoarseNodes->coarseGIDs[currentIndex] = myCoarseGID;
948 coordinatesMap->getRemoteIndexList(ghostedCoarseNodes->GIDs(),
949 ghostedCoarseNodes->PIDs(),
950 ghostedCoarseNodes->LIDs());
961 ghostGIDs.resize(lNumGhostNodes);
964 GO startingGID = minGlobalIndex;
965 Array<GO> startingIndices(3);
968 if(ghostInterface[4] && (myOffset[2] == 0)) {
969 if(gIndices[2] + coarseRate[2] > gFineNodesPerDir[2]) {
970 startingGID -= endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
972 startingGID -= coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
975 startingGID -= myOffset[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
977 if(ghostInterface[2] && (myOffset[1] == 0)) {
978 if(gIndices[1] + coarseRate[1] > gFineNodesPerDir[1]) {
979 startingGID -= endRate[1]*gFineNodesPerDir[0];
981 startingGID -= coarseRate[1]*gFineNodesPerDir[0];
984 startingGID -= myOffset[1]*gFineNodesPerDir[0];
986 if(ghostInterface[0] && (myOffset[0] == 0)) {
987 if(gIndices[0] + coarseRate[0] > gFineNodesPerDir[0]) {
988 startingGID -= endRate[0];
990 startingGID -= coarseRate[0];
993 startingGID -= myOffset[0];
998 startingIndices[2] = startingGID / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
999 tmp = startingGID % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
1000 startingIndices[1] = tmp / gFineNodesPerDir[0];
1001 startingIndices[0] = tmp % gFineNodesPerDir[0];
1004 GO ghostOffset[3] = {0, 0, 0};
1005 LO lengthZero = lCoarseNodesPerDir[0];
1006 LO lengthOne = lCoarseNodesPerDir[1];
1007 LO lengthTwo = lCoarseNodesPerDir[2];
1008 if(ghostInterface[0]) {++lengthZero;}
1009 if(ghostInterface[1]) {++lengthZero;}
1010 if(ghostInterface[2]) {++lengthOne;}
1011 if(ghostInterface[3]) {++lengthOne;}
1012 if(ghostInterface[4]) {++lengthTwo;}
1013 if(ghostInterface[5]) {++lengthTwo;}
1016 if(ghostInterface[4]) {
1017 ghostOffset[2] = startingGID;
1018 for(LO j = 0; j < lengthOne; ++j) {
1019 if( (j == lengthOne-1)
1020 && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1021 ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
1023 ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1025 for(LO k = 0; k < lengthZero; ++k) {
1026 if( (k == lengthZero-1)
1027 && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1028 ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1030 ghostOffset[0] = k*coarseRate[0];
1035 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1043 for(LO i = 0; i < lengthTwo; ++i) {
1046 if( !((i == lengthTwo-1) && ghostInterface[5]) && !((i == 0) && ghostInterface[4]) ) {
1049 if( (i == lengthTwo-1) && !ghostInterface[5] ) {
1050 ghostOffset[2] = startingGID + ((i-1)*coarseRate[2] + endRate[2])
1051 *gFineNodesPerDir[1]*gFineNodesPerDir[0];
1053 ghostOffset[2] = startingGID + i*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1055 for(LO j = 0; j < lengthOne; ++j) {
1056 if( (j == 0) && ghostInterface[2] ) {
1057 for(LO k = 0; k < lengthZero; ++k) {
1058 if( (k == lengthZero-1)
1059 && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1061 ghostOffset[0] = endRate[0];
1063 ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1066 ghostOffset[0] = k*coarseRate[0];
1069 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[0];
1072 }
else if( (j == lengthOne-1) && ghostInterface[3] ) {
1075 if( (j == lengthOne-1)
1076 && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1077 ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
1079 ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1081 for(LO k = 0; k < lengthZero; ++k) {
1082 if( (k == lengthZero-1)
1083 && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1084 ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1086 ghostOffset[0] = k*coarseRate[0];
1088 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1093 if( (j == lengthOne-1)
1094 && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1095 ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
1097 ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1101 if(ghostInterface[0]) {
1102 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1];
1105 if(ghostInterface[1]) {
1106 if( (startingIndices[0] + (lengthZero-1)*coarseRate[0]) > gFineNodesPerDir[0] - 1 ) {
1107 if(lengthZero > 2) {
1108 ghostOffset[0] = (lengthZero-2)*coarseRate[0] + endRate[0];
1110 ghostOffset[0] = endRate[0];
1113 ghostOffset[0] = (lengthZero-1)*coarseRate[0];
1115 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1124 if(ghostInterface[5]) {
1125 if( startingIndices[2] + (lengthTwo-1)*coarseRate[2] + 1 > gFineNodesPerDir[2] ) {
1126 ghostOffset[2] = startingGID
1127 + ((lengthTwo-2)*coarseRate[2] + endRate[2])*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1129 ghostOffset[2] = startingGID
1130 + (lengthTwo-1)*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1132 for(LO j = 0; j < lengthOne; ++j) {
1133 if( (j == lengthOne-1)
1134 && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1135 ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
1137 ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1139 for(LO k = 0; k < lengthZero; ++k) {
1140 if( (k == lengthZero-1)
1141 && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1142 ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1144 ghostOffset[0] = k*coarseRate[0];
1146 ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1159 Array<GO> gCoarseNodesGIDs(lNumCoarseNodes);
1160 LO currentNode, offset2, offset1, offset0;
1162 for(LO ind2 = 0; ind2 < lCoarseNodesPerDir[2]; ++ind2) {
1163 if(myOffset[2] == 0) {
1164 offset2 = startingIndices[2] + myOffset[2];
1166 if(startingIndices[2] + endRate[2] == gFineNodesPerDir[2] - 1) {
1167 offset2 = startingIndices[2] + endRate[2];
1169 offset2 = startingIndices[2] + coarseRate[2];
1172 if(offset2 + ind2*coarseRate[2] > gFineNodesPerDir[2] - 1) {
1173 offset2 += (ind2 - 1)*coarseRate[2] + endRate[2];
1175 offset2 += ind2*coarseRate[2];
1177 offset2 = offset2*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1179 for(LO ind1 = 0; ind1 < lCoarseNodesPerDir[1]; ++ind1) {
1180 if(myOffset[1] == 0) {
1181 offset1 = startingIndices[1] + myOffset[1];
1183 if(startingIndices[1] + endRate[1] == gFineNodesPerDir[1] - 1) {
1184 offset1 = startingIndices[1] + endRate[1];
1186 offset1 = startingIndices[1] + coarseRate[1];
1189 if(offset1 + ind1*coarseRate[1] > gFineNodesPerDir[1] - 1) {
1190 offset1 += (ind1 - 1)*coarseRate[1] + endRate[1];
1192 offset1 += ind1*coarseRate[1];
1194 offset1 = offset1*gFineNodesPerDir[0];
1195 for(LO ind0 = 0; ind0 < lCoarseNodesPerDir[0]; ++ind0) {
1196 offset0 = startingIndices[0];
1197 if(myOffset[0] == 0) {
1198 offset0 += ind0*coarseRate[0];
1200 offset0 += (ind0 + 1)*coarseRate[0];
1202 if(offset0 > gFineNodesPerDir[0] - 1) {offset0 += endRate[0] - coarseRate[0];}
1204 currentNode = ind2*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]
1205 + ind1*lCoarseNodesPerDir[0]
1207 gCoarseNodesGIDs[currentNode] = offset2 + offset1 + offset0;
1214 colGIDs.resize(BlkSize*(lNumCoarseNodes+lNumGhostNodes));
1215 coarseNodesGIDs.resize(lNumCoarseNodes);
1216 for(LO i = 0; i < numDimensions; ++i) {coarseNodes[i].resize(lNumCoarseNodes);}
1217 GO fineNodesPerCoarseSlab = coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1218 GO fineNodesEndCoarseSlab = endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1219 GO fineNodesPerCoarsePlane = coarseRate[1]*gFineNodesPerDir[0];
1220 GO fineNodesEndCoarsePlane = endRate[1]*gFineNodesPerDir[0];
1221 GO coarseNodesPerCoarseLayer = gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0];
1222 GO gCoarseNodeOnCoarseGridGID;
1224 Array<int> ghostPIDs (lNumGhostNodes);
1225 Array<LO> ghostLIDs (lNumGhostNodes);
1226 Array<LO> ghostPermut(lNumGhostNodes);
1227 for(LO k = 0; k < lNumGhostNodes; ++k) {ghostPermut[k] = k;}
1228 coordinatesMap->getRemoteIndexList(ghostGIDs, ghostPIDs, ghostLIDs);
1229 sh_sort_permute(ghostPIDs.begin(),ghostPIDs.end(), ghostPermut.begin(),ghostPermut.end());
1232 GO tmpInds[3], tmpVars[2];
1240 LO firstCoarseNodeInds[3], currentCoarseNode;
1241 for(LO dim = 0; dim < 3; ++dim) {
1242 if(myOffset[dim] == 0) {
1243 firstCoarseNodeInds[dim] = 0;
1245 firstCoarseNodeInds[dim] = coarseRate[dim] - myOffset[dim];
1248 Array<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::magnitudeType> > fineNodes(numDimensions);
1249 for(LO dim = 0; dim < numDimensions; ++dim) {fineNodes[dim] = coordinates->getData(dim);}
1250 for(LO k = 0; k < lCoarseNodesPerDir[2]; ++k) {
1251 for(LO j = 0; j < lCoarseNodesPerDir[1]; ++j) {
1252 for(LO i = 0; i < lCoarseNodesPerDir[0]; ++i) {
1253 col = k*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0] + j*lCoarseNodesPerDir[0] + i;
1256 currentCoarseNode = 0;
1257 if(firstCoarseNodeInds[0] + i*coarseRate[0] > lFineNodesPerDir[0] - 1) {
1258 currentCoarseNode += firstCoarseNodeInds[0] + (i-1)*coarseRate[0] + endRate[0];
1260 currentCoarseNode += firstCoarseNodeInds[0] + i*coarseRate[0];
1262 if(firstCoarseNodeInds[1] + j*coarseRate[1] > lFineNodesPerDir[1] - 1) {
1263 currentCoarseNode += (firstCoarseNodeInds[1] + (j-1)*coarseRate[1] + endRate[1])
1264 *lFineNodesPerDir[0];
1266 currentCoarseNode += (firstCoarseNodeInds[1] + j*coarseRate[1])*lFineNodesPerDir[0];
1268 if(firstCoarseNodeInds[2] + k*coarseRate[2] > lFineNodesPerDir[2] - 1) {
1269 currentCoarseNode += (firstCoarseNodeInds[2] + (k-1)*coarseRate[2] + endRate[2])
1270 *lFineNodesPerDir[1]*lFineNodesPerDir[0];
1272 currentCoarseNode += (firstCoarseNodeInds[2] + k*coarseRate[2])
1273 *lFineNodesPerDir[1]*lFineNodesPerDir[0];
1276 for(LO dim = 0; dim < numDimensions; ++dim) {
1277 coarseNodes[dim][col] = fineNodes[dim][currentCoarseNode];
1280 if((endRate[2] != coarseRate[2])
1281 && (gCoarseNodesGIDs[col] > (gCoarseNodesPerDir[2] - 2)*fineNodesPerCoarseSlab
1282 + fineNodesEndCoarseSlab - 1)) {
1283 tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab + 1;
1284 tmpVars[0] = gCoarseNodesGIDs[col] - (tmpInds[2] - 1)*fineNodesPerCoarseSlab
1285 - fineNodesEndCoarseSlab;
1287 tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab;
1288 tmpVars[0] = gCoarseNodesGIDs[col] % fineNodesPerCoarseSlab;
1290 if((endRate[1] != coarseRate[1])
1291 && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane
1292 + fineNodesEndCoarsePlane - 1)) {
1293 tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
1294 tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane
1295 - fineNodesEndCoarsePlane;
1297 tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
1298 tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
1300 if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
1301 tmpInds[0] = gCoarseNodesPerDir[0] - 1;
1303 tmpInds[0] = tmpVars[1] / coarseRate[0];
1305 gInd[2] = col / (lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]);
1306 tmp = col % (lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]);
1307 gInd[1] = tmp / lCoarseNodesPerDir[0];
1308 gInd[0] = tmp % lCoarseNodesPerDir[0];
1309 lCol = gInd[2]*(lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0])
1310 + gInd[1]*lCoarseNodesPerDir[0] + gInd[0];
1311 gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer
1312 + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
1313 coarseNodesGIDs[lCol] = gCoarseNodeOnCoarseGridGID;
1314 for(LO dof = 0; dof < BlkSize; ++dof) {
1315 colGIDs[BlkSize*lCol + dof] = BlkSize*gCoarseNodeOnCoarseGridGID + dof;
1322 for(col = lNumCoarseNodes; col < lNumCoarseNodes + lNumGhostNodes; ++col) {
1323 if((endRate[2] != coarseRate[2])
1324 && (ghostGIDs[ghostPermut[col - lNumCoarseNodes]] > (gCoarseNodesPerDir[2] - 2)
1325 *fineNodesPerCoarseSlab + fineNodesEndCoarseSlab - 1)) {
1326 tmpInds[2] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab + 1;
1327 tmpVars[0] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]]
1328 - (tmpInds[2] - 1)*fineNodesPerCoarseSlab - fineNodesEndCoarseSlab;
1330 tmpInds[2] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab;
1331 tmpVars[0] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] % fineNodesPerCoarseSlab;
1333 if((endRate[1] != coarseRate[1])
1334 && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane
1335 + fineNodesEndCoarsePlane - 1)) {
1336 tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
1337 tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane
1338 - fineNodesEndCoarsePlane;
1340 tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
1341 tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
1343 if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
1344 tmpInds[0] = gCoarseNodesPerDir[0] - 1;
1346 tmpInds[0] = tmpVars[1] / coarseRate[0];
1348 gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer
1349 + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
1350 for(LO dof = 0; dof < BlkSize; ++dof) {
1351 colGIDs[BlkSize*col + dof] = BlkSize*gCoarseNodeOnCoarseGridGID + dof;
1361 const Array<LO> ,
const LO BlkSize,
const Array<LO> elemInds,
1362 const Array<LO> ,
const LO numDimensions,
1363 const Array<LO> lFineNodesPerDir,
const Array<GO> ,
1364 const Array<GO> ,
const Array<LO> ,
1365 const Array<bool> ghostInterface,
const Array<int> elementFlags,
1366 const std::string stencilType,
const std::string ,
1367 const Array<LO> elementNodesPerDir,
const LO numNodesInElement,
1369 Teuchos::SerialDenseMatrix<LO,SC>& Pi, Teuchos::SerialDenseMatrix<LO,SC>& Pf,
1370 Teuchos::SerialDenseMatrix<LO,SC>& Pe, Array<LO>& dofType,
1371 Array<LO>& lDofInd)
const {
1382 LO countInterior=0, countFace=0, countEdge=0, countCorner =0;
1383 Array<LO> collapseDir(numNodesInElement*BlkSize);
1384 for(LO ke = 0; ke < elementNodesPerDir[2]; ++ke) {
1385 for(LO je = 0; je < elementNodesPerDir[1]; ++je) {
1386 for(LO ie = 0; ie < elementNodesPerDir[0]; ++ie) {
1388 if((ke == 0 || ke == elementNodesPerDir[2]-1)
1389 && (je == 0 || je == elementNodesPerDir[1]-1)
1390 && (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1391 for(LO dof = 0; dof < BlkSize; ++dof) {
1392 dofType[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1393 + je*elementNodesPerDir[0] + ie) + dof] = 0;
1394 lDofInd[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1395 + je*elementNodesPerDir[0] + ie) + dof] = BlkSize*countCorner + dof;
1400 }
else if (((ke == 0 || ke == elementNodesPerDir[2]-1)
1401 && (je == 0 || je == elementNodesPerDir[1]-1))
1402 || ((ke == 0 || ke == elementNodesPerDir[2]-1)
1403 && (ie == 0 || ie == elementNodesPerDir[0]-1))
1404 || ((je == 0 || je == elementNodesPerDir[1]-1)
1405 && (ie == 0 || ie == elementNodesPerDir[0]-1))) {
1406 for(LO dof = 0; dof < BlkSize; ++dof) {
1407 dofType[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1408 + je*elementNodesPerDir[0] + ie) + dof] = 1;
1409 lDofInd[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1410 + je*elementNodesPerDir[0] + ie) + dof] = BlkSize*countEdge + dof;
1411 if((ke == 0 || ke == elementNodesPerDir[2]-1)
1412 && (je == 0 || je == elementNodesPerDir[1]-1)) {
1413 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1414 + je*elementNodesPerDir[0] + ie) + dof] = 0;
1415 }
else if((ke == 0 || ke == elementNodesPerDir[2]-1)
1416 && (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1417 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1418 + je*elementNodesPerDir[0] + ie) + dof] = 1;
1419 }
else if((je == 0 || je == elementNodesPerDir[1]-1
1420 ) && (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1421 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1422 + je*elementNodesPerDir[0] + ie) + dof] = 2;
1428 }
else if ((ke == 0 || ke == elementNodesPerDir[2]-1)
1429 || (je == 0 || je == elementNodesPerDir[1]-1)
1430 || (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1431 for(LO dof = 0; dof < BlkSize; ++dof) {
1432 dofType[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1433 + je*elementNodesPerDir[0] + ie) + dof] = 2;
1434 lDofInd[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1435 + je*elementNodesPerDir[0] + ie) + dof] = BlkSize*countFace + dof;
1436 if(ke == 0 || ke == elementNodesPerDir[2]-1) {
1437 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1438 + je*elementNodesPerDir[0] + ie) + dof] = 2;
1439 }
else if(je == 0 || je == elementNodesPerDir[1]-1) {
1440 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1441 + je*elementNodesPerDir[0] + ie) + dof] = 1;
1442 }
else if(ie == 0 || ie == elementNodesPerDir[0]-1) {
1443 collapseDir[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1444 + je*elementNodesPerDir[0] + ie) + dof] = 0;
1451 for(LO dof = 0; dof < BlkSize; ++dof) {
1452 dofType[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1453 + je*elementNodesPerDir[0] + ie) + dof] = 3;
1454 lDofInd[BlkSize*(ke*elementNodesPerDir[1]*elementNodesPerDir[0]
1455 + je*elementNodesPerDir[0] + ie) + dof] = BlkSize*countInterior +dof;
1463 LO numInteriorNodes = 0, numFaceNodes = 0, numEdgeNodes = 0, numCornerNodes = 8;
1464 numInteriorNodes = (elementNodesPerDir[0] - 2)*(elementNodesPerDir[1] - 2)
1465 *(elementNodesPerDir[2] - 2);
1466 numFaceNodes = 2*(elementNodesPerDir[0] - 2)*(elementNodesPerDir[1] - 2)
1467 + 2*(elementNodesPerDir[0] - 2)*(elementNodesPerDir[2] - 2)
1468 + 2*(elementNodesPerDir[1] - 2)*(elementNodesPerDir[2] - 2);
1469 numEdgeNodes = 4*(elementNodesPerDir[0] - 2) + 4*(elementNodesPerDir[1] - 2)
1470 + 4*(elementNodesPerDir[2] - 2);
1472 Teuchos::SerialDenseMatrix<LO,SC> Aii(BlkSize*numInteriorNodes, BlkSize*numInteriorNodes);
1473 Teuchos::SerialDenseMatrix<LO,SC> Aff(BlkSize*numFaceNodes, BlkSize*numFaceNodes);
1474 Teuchos::SerialDenseMatrix<LO,SC> Aee(BlkSize*numEdgeNodes, BlkSize*numEdgeNodes);
1476 Teuchos::SerialDenseMatrix<LO,SC> Aif(BlkSize*numInteriorNodes, BlkSize*numFaceNodes);
1477 Teuchos::SerialDenseMatrix<LO,SC> Aie(BlkSize*numInteriorNodes, BlkSize*numEdgeNodes);
1478 Teuchos::SerialDenseMatrix<LO,SC> Afe(BlkSize*numFaceNodes, BlkSize*numEdgeNodes);
1480 Teuchos::SerialDenseMatrix<LO,SC> Aic(BlkSize*numInteriorNodes, BlkSize*numCornerNodes);
1481 Teuchos::SerialDenseMatrix<LO,SC> Afc(BlkSize*numFaceNodes, BlkSize*numCornerNodes);
1482 Teuchos::SerialDenseMatrix<LO,SC> Aec(BlkSize*numEdgeNodes, BlkSize*numCornerNodes);
1483 Pi.shape( BlkSize*numInteriorNodes, BlkSize*numCornerNodes);
1484 Pf.shape( BlkSize*numFaceNodes, BlkSize*numCornerNodes);
1485 Pe.shape( BlkSize*numEdgeNodes, BlkSize*numCornerNodes);
1487 ArrayView<const LO> rowIndices;
1488 ArrayView<const SC> rowValues;
1489 LO idof, iInd, jInd;
1490 int iType = 0, jType = 0;
1491 int orientation = -1;
1492 int collapseFlags[3] = {};
1493 Array<SC> stencil((std::pow(3,numDimensions))*BlkSize);
1497 for(LO ke = 0; ke < elementNodesPerDir[2]; ++ke) {
1498 for(LO je = 0; je < elementNodesPerDir[1]; ++je) {
1499 for(LO ie = 0; ie < elementNodesPerDir[0]; ++ie) {
1500 collapseFlags[0] = 0; collapseFlags[1] = 0; collapseFlags[2] = 0;
1501 if((elementFlags[0] == 1 || elementFlags[0] == 3) && ie == 0) {
1502 collapseFlags[0] += 1;
1504 if((elementFlags[0] == 2 || elementFlags[0] == 3) && ie == elementNodesPerDir[0] - 1) {
1505 collapseFlags[0] += 2;
1507 if((elementFlags[1] == 1 || elementFlags[1] == 3) && je == 0) {
1508 collapseFlags[1] += 1;
1510 if((elementFlags[1] == 2 || elementFlags[1] == 3) && je == elementNodesPerDir[1] - 1) {
1511 collapseFlags[1] += 2;
1513 if((elementFlags[2] == 1 || elementFlags[2] == 3) && ke == 0) {
1514 collapseFlags[2] += 1;
1516 if((elementFlags[2] == 2 || elementFlags[2] == 3) && ke == elementNodesPerDir[2] - 1) {
1517 collapseFlags[2] += 2;
1521 GetNodeInfo(ie, je, ke, elementNodesPerDir, &iType, iInd, &orientation);
1522 for(LO dof0 = 0; dof0 < BlkSize; ++dof0) {
1524 idof = ((elemInds[2]*coarseRate[2] + ke)*lFineNodesPerDir[1]*lFineNodesPerDir[0]
1525 + (elemInds[1]*coarseRate[1] + je)*lFineNodesPerDir[0]
1526 + elemInds[0]*coarseRate[0] + ie)*BlkSize + dof0;
1527 Aghost->getLocalRowView(idof, rowIndices, rowValues);
1528 FormatStencil(BlkSize, ghostInterface, ie, je, ke, rowValues,
1529 elementNodesPerDir, collapseFlags, stencilType, stencil);
1532 for(LO interactingNode = 0; interactingNode < 27; ++interactingNode) {
1534 ko = ke + (interactingNode / 9 - 1);
1536 LO tmp = interactingNode % 9;
1537 jo = je + (tmp / 3 - 1);
1538 io = ie + (tmp % 3 - 1);
1540 GetNodeInfo(io, jo, ko, elementNodesPerDir, &jType, jInd, &dummy);
1542 if((io > -1 && io < elementNodesPerDir[0]) &&
1543 (jo > -1 && jo < elementNodesPerDir[1]) &&
1544 (ko > -1 && ko < elementNodesPerDir[2])) {
1545 for(LO dof1 = 0; dof1 < BlkSize; ++dof1) {
1547 Aii(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1548 stencil[interactingNode*BlkSize + dof1];
1549 }
else if(jType == 2) {
1550 Aif(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1551 stencil[interactingNode*BlkSize + dof1];
1552 }
else if(jType == 1) {
1553 Aie(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1554 stencil[interactingNode*BlkSize + dof1];
1555 }
else if(jType == 0) {
1556 Aic(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1557 -stencil[interactingNode*BlkSize + dof1];
1562 }
else if(iType == 2) {
1563 CollapseStencil(2, orientation, collapseFlags, stencil);
1564 for(LO interactingNode = 0; interactingNode < 27; ++interactingNode) {
1566 ko = ke + (interactingNode / 9 - 1);
1568 LO tmp = interactingNode % 9;
1569 jo = je + (tmp / 3 - 1);
1570 io = ie + (tmp % 3 - 1);
1572 GetNodeInfo(io, jo, ko, elementNodesPerDir, &jType, jInd, &dummy);
1574 if((io > -1 && io < elementNodesPerDir[0]) &&
1575 (jo > -1 && jo < elementNodesPerDir[1]) &&
1576 (ko > -1 && ko < elementNodesPerDir[2])) {
1577 for(LO dof1 = 0; dof1 < BlkSize; ++dof1) {
1579 Aff(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1580 stencil[interactingNode*BlkSize + dof1];
1581 }
else if(jType == 1) {
1582 Afe(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1583 stencil[interactingNode*BlkSize + dof1];
1584 }
else if(jType == 0) {
1585 Afc(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1586 -stencil[interactingNode*BlkSize + dof1];
1591 }
else if(iType == 1) {
1592 CollapseStencil(1, orientation, collapseFlags, stencil);
1593 for(LO interactingNode = 0; interactingNode < 27; ++interactingNode) {
1595 ko = ke + (interactingNode / 9 - 1);
1597 LO tmp = interactingNode % 9;
1598 jo = je + (tmp / 3 - 1);
1599 io = ie + (tmp % 3 - 1);
1601 GetNodeInfo(io, jo, ko, elementNodesPerDir, &jType, jInd, &dummy);
1603 if((io > -1 && io < elementNodesPerDir[0]) &&
1604 (jo > -1 && jo < elementNodesPerDir[1]) &&
1605 (ko > -1 && ko < elementNodesPerDir[2])) {
1606 for(LO dof1 = 0; dof1 < BlkSize; ++dof1) {
1608 Aee(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1609 stencil[interactingNode*BlkSize + dof1];
1610 }
else if(jType == 0) {
1611 Aec(iInd*BlkSize + dof0, jInd*BlkSize + dof1) =
1612 -stencil[interactingNode*BlkSize + dof1];
1633 Teuchos::SerialDenseSolver<LO,SC> problem;
1634 problem.setMatrix(Teuchos::rcp(&Aee,
false));
1635 problem.setVectors(Teuchos::rcp(&Pe,
false), Teuchos::rcp(&Aec,
false));
1636 problem.factorWithEquilibration(
true);
1638 problem.unequilibrateLHS();
1644 Afc.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -1.0, Afe, Pe, 1.0);
1646 Teuchos::SerialDenseSolver<LO,SC> problem;
1647 problem.setMatrix(Teuchos::rcp(&Aff,
false));
1648 problem.setVectors(Teuchos::rcp(&Pf,
false), Teuchos::rcp(&Afc,
false));
1649 problem.factorWithEquilibration(
true);
1651 problem.unequilibrateLHS();
1657 Aic.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -1.0, Aie, Pe, 1.0);
1659 Aic.multiply(Teuchos::NO_TRANS, Teuchos::NO_TRANS, -1.0, Aif, Pf, 1.0);
1661 Teuchos::SerialDenseSolver<LO,SC> problem;
1662 problem.setMatrix(Teuchos::rcp(&Aii,
false));
1663 problem.setVectors(Teuchos::rcp(&Pi,
false), Teuchos::rcp(&Aic,
false));
1664 problem.factorWithEquilibration(
true);
1666 problem.unequilibrateLHS();