46#ifndef MUELU_BLACKBOXPFACTORY_DEF_HPP
47#define MUELU_BLACKBOXPFACTORY_DEF_HPP
54#include <Teuchos_SerialDenseMatrix.hpp>
55#include <Teuchos_SerialDenseVector.hpp>
56#include <Teuchos_SerialDenseSolver.hpp>
58#include <Xpetra_CrsMatrixUtils.hpp>
59#include <Xpetra_CrsMatrixWrap.hpp>
60#include <Xpetra_ImportFactory.hpp>
61#include <Xpetra_Matrix.hpp>
62#include <Xpetra_MapFactory.hpp>
63#include <Xpetra_MultiVectorFactory.hpp>
64#include <Xpetra_VectorFactory.hpp>
66#include <Xpetra_IO.hpp>
76 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 RCP<ParameterList> validParamList = rcp(
new ParameterList());
85 validParamList->set<std::string > (
"Coarsen",
"{3}",
"Coarsening rate in all spatial dimensions");
86 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
87 validParamList->set<RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
"Generating factory of the nullspace");
88 validParamList->set<RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for coorindates");
89 validParamList->set<RCP<const FactoryBase> >(
"gNodesPerDim", Teuchos::null,
"Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
90 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
"Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
91 validParamList->set<std::string> (
"stencil type",
"full",
"You can use two type of stencils: full and reduced, that correspond to 27 and 7 points stencils respectively in 3D.");
92 validParamList->set<std::string> (
"block strategy",
"coupled",
"The strategy used to handle systems of PDEs can be: coupled or uncoupled.");
94 return validParamList;
97 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
101 Input(fineLevel,
"A");
102 Input(fineLevel,
"Nullspace");
103 Input(fineLevel,
"Coordinates");
111 "gNodesPerDim was not provided by the user on level0!");
114 Input(fineLevel,
"gNodesPerDim");
124 "lNodesPerDim was not provided by the user on level0!");
127 Input(fineLevel,
"lNodesPerDim");
131 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
133 Level& coarseLevel)
const{
134 return BuildP(fineLevel, coarseLevel);
138 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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);
757 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
1358 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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();
1671 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1674 Array<SC>& stencil)
const {
1678 if(orientation == 0) {
1679 for(LO i = 0; i < 9; ++i) {
1680 stencil[3*i + 1] = stencil[3*i] + stencil[3*i + 1] + stencil[3*i + 2];
1682 stencil[3*i + 2] = 0;
1684 }
else if(orientation == 1) {
1685 for(LO i = 0; i < 3; ++i) {
1686 stencil[9*i + 3] = stencil[9*i + 0] + stencil[9*i + 3] + stencil[9*i + 6];
1687 stencil[9*i + 0] = 0;
1688 stencil[9*i + 6] = 0;
1689 stencil[9*i + 4] = stencil[9*i + 1] + stencil[9*i + 4] + stencil[9*i + 7];
1690 stencil[9*i + 1] = 0;
1691 stencil[9*i + 7] = 0;
1692 stencil[9*i + 5] = stencil[9*i + 2] + stencil[9*i + 5] + stencil[9*i + 8];
1693 stencil[9*i + 2] = 0;
1694 stencil[9*i + 8] = 0;
1696 }
else if(orientation == 2) {
1697 for(LO i = 0; i < 9; ++i) {
1698 stencil[i + 9] = stencil[i] + stencil[i + 9] + stencil[i + 18];
1700 stencil[i + 18] = 0;
1703 }
else if(type == 1) {
1704 SC tmp1 = 0, tmp2 = 0, tmp3 = 0;
1705 if(orientation == 0) {
1706 for(LO i = 0; i < 9; ++i) {
1707 tmp1 += stencil[0 + 3*i];
1708 stencil[0 + 3*i] = 0;
1709 tmp2 += stencil[1 + 3*i];
1710 stencil[1 + 3*i] = 0;
1711 tmp3 += stencil[2 + 3*i];
1712 stencil[2 + 3*i] = 0;
1717 }
else if(orientation == 1) {
1718 for(LO i = 0; i < 3; ++i) {
1719 tmp1 += stencil[0 + i] + stencil[9 + i] + stencil[18 + i];
1720 stencil[ 0 + i] = 0;
1721 stencil[ 9 + i] = 0;
1722 stencil[18 + i] = 0;
1723 tmp2 += stencil[3 + i] + stencil[12 + i] + stencil[21 + i];
1724 stencil[ 3 + i] = 0;
1725 stencil[12 + i] = 0;
1726 stencil[21 + i] = 0;
1727 tmp3 += stencil[6 + i] + stencil[15 + i] + stencil[24 + i];
1728 stencil[ 6 + i] = 0;
1729 stencil[15 + i] = 0;
1730 stencil[24 + i] = 0;
1735 }
else if(orientation == 2) {
1736 for(LO i = 0; i < 9; ++i) {
1739 tmp2 += stencil[i + 9];
1741 tmp3 += stencil[i + 18];
1742 stencil[i + 18] = 0;
1751 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1753 FormatStencil(
const LO BlkSize,
const Array<bool> ,
const LO ,
const LO ,
1754 const LO ,
const ArrayView<const SC> rowValues,
const Array<LO> ,
1755 const int collapseFlags[3],
const std::string stencilType, Array<SC>& stencil)
1758 if(stencilType ==
"reduced") {
1759 Array<int> fullStencilInds(7);
1760 fullStencilInds[0] = 4; fullStencilInds[1] = 10; fullStencilInds[2] = 12;
1761 fullStencilInds[3] = 13; fullStencilInds[4] = 14; fullStencilInds[5] = 16;
1762 fullStencilInds[6] = 22;
1766 int stencilSize = rowValues.size();
1767 if(collapseFlags[0] + collapseFlags[1] + collapseFlags[2] > 0) {
1768 if(stencilSize == 1) {
1772 mask[0] = 1; mask[1] = 1; mask[2] = 1;
1773 mask[4] = 1; mask[5] = 1; mask[6] = 1;
1777 if(collapseFlags[2] == 1 || collapseFlags[2] == 3) {
1780 if(collapseFlags[2] == 2 || collapseFlags[2] == 3) {
1783 if(collapseFlags[1] == 1 || collapseFlags[1] == 3) {
1786 if(collapseFlags[1] == 2 || collapseFlags[1] == 3) {
1789 if(collapseFlags[0] == 1 || collapseFlags[0] == 3) {
1792 if(collapseFlags[0] == 2 || collapseFlags[0] == 3) {
1799 for(
int ind = 0; ind < 7; ++ind) {
1800 if(mask[ind] == 1) {
1801 for(
int dof = 0; dof < BlkSize; ++dof) {
1802 stencil[BlkSize*fullStencilInds[ind] + dof] = 0.0;
1806 for(
int dof = 0; dof < BlkSize; ++dof) {
1807 stencil[BlkSize*fullStencilInds[ind] + dof] = rowValues[BlkSize*(ind - offset) + dof];
1811 }
else if(stencilType ==
"full") {
1813 Array<int> mask(27);
1814 if(collapseFlags[2] == 1 || collapseFlags[2] == 3) {
1815 for(
int ind = 0; ind < 9; ++ind) {
1819 if(collapseFlags[2] == 2 || collapseFlags[2] == 3) {
1820 for(
int ind = 0; ind < 9; ++ind) {
1824 if(collapseFlags[1] == 1 || collapseFlags[1] == 3) {
1825 for(
int ind1 = 0; ind1 < 3; ++ind1) {
1826 for(
int ind2 = 0; ind2 < 3; ++ind2) {
1827 mask[ind1*9 + ind2] = 1;
1831 if(collapseFlags[1] == 2 || collapseFlags[1] == 3) {
1832 for(
int ind1 = 0; ind1 < 3; ++ind1) {
1833 for(
int ind2 = 0; ind2 < 3; ++ind2) {
1834 mask[ind1*9 + ind2 + 6] = 1;
1838 if(collapseFlags[0] == 1 || collapseFlags[0] == 3) {
1839 for(
int ind = 0; ind < 9; ++ind) {
1843 if(collapseFlags[0] == 2 || collapseFlags[0] == 3) {
1844 for(
int ind = 0; ind < 9; ++ind) {
1845 mask[3*ind + 2] = 1;
1851 for(LO ind = 0; ind < 27; ++ind) {
1852 if(mask[ind] == 0) {
1853 for(
int dof = 0; dof < BlkSize; ++dof) {
1854 stencil[BlkSize*ind + dof] = 0.0;
1858 for(
int dof = 0; dof < BlkSize; ++dof) {
1859 stencil[BlkSize*ind + dof] = rowValues[BlkSize*(ind - offset) + dof];
1866 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1868 const LO ie,
const LO je,
const LO ke,
1869 const Array<LO> elementNodesPerDir,
1870 int* type, LO& ind,
int* orientation)
const {
1873 *type = -1, ind = 0, *orientation = -1;
1874 if((ke == 0 || ke == elementNodesPerDir[2]-1)
1875 && (je == 0 || je == elementNodesPerDir[1]-1)
1876 && (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1879 ind = 4*ke / (elementNodesPerDir[2]-1) + 2*je / (elementNodesPerDir[1]-1)
1880 + ie / (elementNodesPerDir[0]-1);
1881 }
else if(((ke == 0 || ke == elementNodesPerDir[2]-1)
1882 && (je == 0 || je == elementNodesPerDir[1]-1))
1883 || ((ke == 0 || ke == elementNodesPerDir[2]-1)
1884 && (ie == 0 || ie == elementNodesPerDir[0]-1))
1885 || ((je == 0 || je == elementNodesPerDir[1]-1)
1886 && (ie == 0 || ie == elementNodesPerDir[0]-1))) {
1889 if(ke > 0) {ind += 2*(elementNodesPerDir[0] - 2 + elementNodesPerDir[1] - 2);}
1890 if(ke == elementNodesPerDir[2] - 1) {ind += 4*(elementNodesPerDir[2] - 2);}
1891 if((ke == 0) || (ke == elementNodesPerDir[2] - 1)) {
1895 }
else if(je == elementNodesPerDir[1] - 1) {
1897 ind += 2*(elementNodesPerDir[1] - 2) + elementNodesPerDir[0] - 2 + ie - 1;
1900 ind += elementNodesPerDir[0] - 2 + 2*(je - 1) + ie / (elementNodesPerDir[0] - 1);
1904 ind += 4*(ke - 1) + 2*(je/(elementNodesPerDir[1] - 1)) + ie / (elementNodesPerDir[0] - 1);
1906 }
else if ((ke == 0 || ke == elementNodesPerDir[2]-1)
1907 || (je == 0 || je == elementNodesPerDir[1]-1)
1908 || (ie == 0 || ie == elementNodesPerDir[0]-1)) {
1913 ind = (je - 1)*(elementNodesPerDir[0] - 2) + ie - 1;
1916 ind += (elementNodesPerDir[1] - 2)*(elementNodesPerDir[0] - 2);
1918 ind += 2*(ke - 1)*(elementNodesPerDir[1] - 2 + elementNodesPerDir[0] - 2);
1919 if(ke == elementNodesPerDir[2]-1) {
1921 ind += (je - 1)*(elementNodesPerDir[0] - 2) + ie - 1;
1926 }
else if(je == elementNodesPerDir[1] - 1) {
1928 ind += 2*(elementNodesPerDir[1] - 2) + elementNodesPerDir[0] - 2 + ie - 1;
1931 ind += elementNodesPerDir[0] - 2 + 2*(je - 1) + ie / (elementNodesPerDir[0]-1);
1938 ind = (ke - 1)*(elementNodesPerDir[1] - 2)*(elementNodesPerDir[0] - 2)
1939 + (je - 1)*(elementNodesPerDir[0] - 2) + ie - 1;
1943 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1945 const typename Teuchos::Array<LocalOrdinal>::iterator& first1,
1946 const typename Teuchos::Array<LocalOrdinal>::iterator& last1,
1947 const typename Teuchos::Array<LocalOrdinal>::iterator& first2,
1948 const typename Teuchos::Array<LocalOrdinal>::iterator& )
const
1950 typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1951 DT n = last1 - first1;
1953 DT z = Teuchos::OrdinalTraits<DT>::zero();
1957 for (DT j = 0; j < max; j++)
1959 for (DT k = j; k >= 0; k-=m)
1961 if (first1[first2[k+m]] >= first1[first2[k]])
1963 std::swap(first2[k+m], first2[k]);
1972#define MUELU_BLACKBOXPFACTORY_SHORT
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void FormatStencil(const LO BlkSize, const Array< bool > ghostInterface, const LO ie, const LO je, const LO ke, const ArrayView< const SC > rowValues, const Array< LO > elementNodesPerDir, const int collapseFlags[3], const std::string stencilType, Array< SC > &stencil) const
void GetNodeInfo(const LO ie, const LO je, const LO ke, const Array< LO > elementNodesPerDir, int *type, LO &ind, int *orientation) const
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
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void ComputeLocalEntries(const RCP< const Matrix > &Aghost, const Array< LO > coarseRate, const Array< LO > endRate, const LO BlkSize, const Array< LO > elemInds, const Array< LO > lCoarseElementsPerDir, const LO numDimensions, const Array< LO > lFineNodesPerDir, const Array< GO > gFineNodesPerDir, const Array< GO > gIndices, const Array< LO > lCoarseNodesPerDir, const Array< bool > ghostInterface, const Array< int > elementFlags, const std::string stencilType, const std::string blockStrategy, const Array< LO > elementNodesPerDir, const LO numNodesInElement, const Array< GO > colGIDs, Teuchos::SerialDenseMatrix< LO, SC > &Pi, Teuchos::SerialDenseMatrix< LO, SC > &Pf, Teuchos::SerialDenseMatrix< LO, SC > &Pe, Array< LO > &dofType, Array< LO > &lDofInd) const
void CollapseStencil(const int type, const int orientation, const int collapseFlags[3], Array< SC > &stencil) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void GetGeometricData(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LO, GO, NO > > &coordinates, const Array< LO > coarseRate, const Array< GO > gFineNodesPerDir, const Array< LO > lFineNodesPerDir, const LO BlkSize, Array< GO > &gIndices, Array< LO > &myOffset, Array< bool > &ghostInterface, Array< LO > &endRate, Array< GO > &gCoarseNodesPerDir, Array< LO > &lCoarseNodesPerDir, Array< LO > &glCoarseNodesPerDir, Array< GO > &ghostGIDs, Array< GO > &coarseNodesGIDs, Array< GO > &colGIDs, GO &gNumCoarseNodes, LO &lNumCoarseNodes, ArrayRCP< Array< typename Teuchos::ScalarTraits< Scalar >::magnitudeType > > coarseNodes, Array< int > &boundaryFlags, RCP< NodesIDs > ghostedCoarseNodes) const
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
static const NoFactory * get()
Namespace for MueLu classes and methods.
@ Runtime0
One-liner description of what is happening.