46#ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
47#define MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
49#include "Xpetra_CrsGraph.hpp"
50#include "Xpetra_CrsMatrixUtils.hpp"
54#include "MueLu_Aggregates.hpp"
55#ifdef HAVE_MUELU_TPETRA
56#include "Xpetra_TpetraCrsMatrix.hpp"
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
66 RCP<ParameterList> validParamList = rcp(
new ParameterList());
68#define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
73 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
74 "Generating factory of the matrix A");
75 validParamList->set<RCP<const FactoryBase> >(
"Aggregates", Teuchos::null,
76 "Aggregates generated by StructuredAggregationFactory used to construct a piece-constant prolongator.");
77 validParamList->set<RCP<const FactoryBase> >(
"prolongatorGraph", Teuchos::null,
78 "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
79 validParamList->set<RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
80 "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
81 validParamList->set<RCP<const FactoryBase> >(
"coarseCoordinatesFineMap", Teuchos::null,
82 "map of the coarse coordinates' GIDs as indexed on the fine mesh.");
83 validParamList->set<RCP<const FactoryBase> >(
"coarseCoordinatesMap", Teuchos::null,
84 "map of the coarse coordinates' GIDs as indexed on the coarse mesh.");
85 validParamList->set<RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
86 "Fine level nullspace used to construct the coarse level nullspace.");
87 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
88 "Number of spacial dimensions in the problem.");
89 validParamList->set<RCP<const FactoryBase> >(
"lCoarseNodesPerDim", Teuchos::null,
90 "Number of nodes per spatial dimension on the coarse grid.");
91 validParamList->set<RCP<const FactoryBase> >(
"structuredInterpolationOrder", Teuchos::null,
92 "Interpolation order for constructing the prolongator.");
93 validParamList->set<
bool> (
"keep coarse coords",
false,
"Flag to keep coordinates for special coarse grid solve");
94 validParamList->set<
bool> (
"interp: remove small entries",
true,
"Remove small interpolation coeficient from prolongator to reduce fill-in on coarse level");
96 return validParamList;
99 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
102 const ParameterList& pL = GetParameterList();
104 Input(fineLevel,
"A");
105 Input(fineLevel,
"Nullspace");
106 Input(fineLevel,
"numDimensions");
107 Input(fineLevel,
"prolongatorGraph");
108 Input(fineLevel,
"lCoarseNodesPerDim");
109 Input(fineLevel,
"structuredInterpolationOrder");
111 if( pL.get<
bool>(
"interp: build coarse coordinates") ||
112 Get<int>(fineLevel,
"structuredInterpolationOrder") == 1) {
113 Input(fineLevel,
"Coordinates");
114 Input(fineLevel,
"coarseCoordinatesFineMap");
115 Input(fineLevel,
"coarseCoordinatesMap");
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 return BuildP(fineLevel, coarseLevel);
125 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
131 RCP<Teuchos::FancyOStream> out;
132 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
133 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
134 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
136 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
139 *out <<
"Starting GeometricInterpolationPFactory::BuildP." << std::endl;
142 const ParameterList& pL = GetParameterList();
143 const bool removeSmallEntries = pL.get<
bool>(
"interp: remove small entries");
144 const bool buildCoarseCoordinates = pL.get<
bool>(
"interp: build coarse coordinates");
145 const int interpolationOrder = Get<int>(fineLevel,
"structuredInterpolationOrder");
146 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
149 RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel,
"A");
150 RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel,
"prolongatorGraph");
151 RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
156 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
158 RCP<const Map> coarseCoordsFineMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesFineMap");
159 RCP<const Map> coarseCoordsMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesMap");
161 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
162 coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::Build(coarseCoordsFineMap,
163 fineCoordinates->getNumVectors());
164 RCP<const Import> coordsImporter = ImportFactory::Build(fineCoordinates->getMap(),
165 coarseCoordsFineMap);
166 coarseCoordinates->doImport(*fineCoordinates, *coordsImporter, Xpetra::INSERT);
167 coarseCoordinates->replaceMap(coarseCoordsMap);
169 Set(coarseLevel,
"Coordinates", coarseCoordinates);
171 if(pL.get<
bool>(
"keep coarse coords")) {
172 coarseLevel.
Set<RCP<realvaluedmultivector_type> >(
"Coordinates2", coarseCoordinates,
NoFactory::get());
176 *out <<
"Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
179 if(interpolationOrder == 0) {
182 BuildConstantP(P, prolongatorGraph, A);
183 }
else if(interpolationOrder == 1) {
186 RCP<const Map> prolongatorColMap = prolongatorGraph->getColMap();
188 const size_t dofsPerNode =
static_cast<size_t>(A->GetFixedBlockSize());
189 const size_t numColIndices = prolongatorColMap->getLocalNumElements();
190 TEUCHOS_TEST_FOR_EXCEPTION((numColIndices % dofsPerNode) != 0,
192 "Something went wrong, the number of columns in the prolongator is not a multiple of dofsPerNode!");
193 const size_t numGhostCoords = numColIndices / dofsPerNode;
194 const GO indexBase = prolongatorColMap->getIndexBase();
195 const GO coordIndexBase = fineCoordinates->getMap()->getIndexBase();
197 ArrayView<const GO> prolongatorColIndices = prolongatorColMap->getLocalElementList();
198 Array<GO> ghostCoordIndices(numGhostCoords);
199 for(
size_t ghostCoordIdx = 0; ghostCoordIdx < numGhostCoords; ++ghostCoordIdx) {
200 ghostCoordIndices[ghostCoordIdx]
201 = (prolongatorColIndices[ghostCoordIdx*dofsPerNode] - indexBase) / dofsPerNode
204 RCP<Map> ghostCoordMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
205 prolongatorColMap->getGlobalNumElements() / dofsPerNode,
208 fineCoordinates->getMap()->getComm());
210 RCP<realvaluedmultivector_type> ghostCoordinates
211 = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(ghostCoordMap,
212 fineCoordinates->getNumVectors());
213 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
215 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
217 BuildLinearP(coarseLevel, A, prolongatorGraph, fineCoordinates,
218 ghostCoordinates, numDimensions, removeSmallEntries, P);
221 *out <<
"The prolongator matrix has been built." << std::endl;
226 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
227 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
228 fineNullspace->getNumVectors());
230 using helpers=Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
231 if(helpers::isTpetraBlockCrs(A)) {
234 Ptrans->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
235 Teuchos::ScalarTraits<SC>::zero());
238 P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
239 Teuchos::ScalarTraits<SC>::zero());
242 Set(coarseLevel,
"Nullspace", coarseNullspace);
245 *out <<
"The coarse nullspace is constructed and set on the coarse level." << std::endl;
247 Set(coarseLevel,
"P", P);
249 *out <<
"GeometricInterpolationPFactory::BuildP has completed." << std::endl;
253 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
255 BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A)
const {
258 RCP<Teuchos::FancyOStream> out;
259 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
260 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
261 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
263 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
266 *out <<
"BuildConstantP" << std::endl;
268 std::vector<size_t> strideInfo(1);
269 strideInfo[0] = A->GetFixedBlockSize();
270 RCP<const StridedMap> stridedDomainMap =
271 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
273 *out <<
"Call prolongator constructor" << std::endl;
275 using helpers=Xpetra::Helpers<Scalar,LocalOrdinal,GlobalOrdinal,Node>;
276 if(helpers::isTpetraBlockCrs(A)) {
277#ifdef HAVE_MUELU_TPETRA
278 SC one = Teuchos::ScalarTraits<SC>::one();
279 SC zero = Teuchos::ScalarTraits<SC>::zero();
280 LO NSDim = A->GetStorageBlockSize();
283 RCP<const Map> BlockMap = prolongatorGraph->getDomainMap();
284 Teuchos::ArrayView<const GO> block_dofs = BlockMap->getLocalElementList();
285 Teuchos::Array<GO> point_dofs(block_dofs.size()*NSDim);
286 for(LO i=0, ct=0; i<block_dofs.size(); i++) {
287 for(LO j=0; j<NSDim; j++) {
288 point_dofs[ct] = block_dofs[i]*NSDim + j;
293 RCP<const Map> PointMap = MapFactory::Build(BlockMap->lib(),
294 BlockMap->getGlobalNumElements() *NSDim,
296 BlockMap->getIndexBase(),
297 BlockMap->getComm());
298 strideInfo[0] = A->GetFixedBlockSize();
299 RCP<const StridedMap> stridedPointMap = StridedMapFactory::Build(PointMap, strideInfo);
301 RCP<Xpetra::CrsMatrix<SC,LO,GO,NO> > P_xpetra = Xpetra::CrsMatrixFactory<SC,LO,GO,NO>::BuildBlock(prolongatorGraph, PointMap, A->getRangeMap(),NSDim);
302 RCP<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> > P_tpetra = rcp_dynamic_cast<Xpetra::TpetraBlockCrsMatrix<SC,LO,GO,NO> >(P_xpetra);
303 if(P_tpetra.is_null())
throw std::runtime_error(
"BuildConstantP Matrix factory did not return a Tpetra::BlockCrsMatrix");
304 RCP<CrsMatrixWrap> P_wrap = rcp(
new CrsMatrixWrap(P_xpetra));
307 Teuchos::Array<LO> temp(1);
308 Teuchos::ArrayView<const LO> indices;
309 Teuchos::Array<Scalar> block(NSDim*NSDim, zero);
310 for(LO i=0; i<NSDim; i++)
311 block[i*NSDim+i] = one;
312 for(LO i=0; i<(int)prolongatorGraph->getLocalNumRows(); i++) {
313 prolongatorGraph->getLocalRowView(i,indices);
314 for(LO j=0; j<(LO)indices.size();j++) {
315 temp[0] = indices[j];
316 P_tpetra->replaceLocalValues(i,temp(),block());
321 if (A->IsView(
"stridedMaps") ==
true) {
322 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedPointMap);
325 P->CreateView(
"stridedMaps", P->getRangeMap(), PointMap);
329 throw std::runtime_error(
"GeometricInteroplationFactory::Build(): BlockCrs requires Tpetra");
335 RCP<ParameterList> dummyList = rcp(
new ParameterList());
336 P = rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
337 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
338 PCrs->setAllToScalar(1.0);
339 PCrs->fillComplete();
342 if (A->IsView(
"stridedMaps") ==
true)
343 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
345 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
350 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
353 RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
354 RCP<realvaluedmultivector_type>& fineCoordinates,
355 RCP<realvaluedmultivector_type>& ghostCoordinates,
356 const int numDimensions,
const bool removeSmallEntries,
357 RCP<Matrix>& P)
const {
360 RCP<Teuchos::FancyOStream> out;
361 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
362 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
363 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
365 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
369 const int numInterpolationPoints = 1 << numDimensions;
370 const int dofsPerNode = A->GetFixedBlockSize()/ A->GetStorageBlockSize();;
372 RCP<ParameterList> dummyList = rcp(
new ParameterList());
373 P = rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
374 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
378 *out <<
"Entering BuildLinearP" << std::endl;
382 const LO numFineNodes = fineCoordinates->getLocalLength();
383 const LO numGhostNodes = ghostCoordinates->getLocalLength();
384 Array<ArrayRCP<const real_type> > fineCoords(3);
385 Array<ArrayRCP<const real_type> > ghostCoords(3);
386 const real_type realZero = Teuchos::as<real_type>(0.0);
387 ArrayRCP<real_type> fineZero(numFineNodes, realZero);
388 ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
389 for(
int dim = 0; dim < 3; ++dim) {
390 if(dim < numDimensions) {
391 fineCoords[dim] = fineCoordinates->getData(dim);
392 ghostCoords[dim] = ghostCoordinates->getData(dim);
394 fineCoords[dim] = fineZero;
395 ghostCoords[dim] = ghostZero;
399 *out <<
"Coordinates extracted from the multivectors!" << std::endl;
402 LO interpolationNodeIdx = 0, rowIdx = 0;
403 ArrayView<const LO> colIndices;
405 Array<Array<real_type> > coords(numInterpolationPoints + 1);
406 Array<real_type> stencil(numInterpolationPoints);
407 for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
408 if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
411 for(LO dof = 0; dof < dofsPerNode; ++dof) {
412 rowIdx = nodeIdx*dofsPerNode + dof;
413 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
414 PCrs->replaceLocalValues(rowIdx, colIndices, values());
420 for(
int dim = 0; dim < 3; ++dim) {
421 coords[0][dim] = fineCoords[dim][nodeIdx];
423 prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
424 for(
int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
425 coords[interpolationIdx + 1].resize(3);
426 interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
427 for(
int dim = 0; dim < 3; ++dim) {
428 coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
431 RCP<Teuchos::TimeMonitor> tm = rcp(
new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(
"Compute Linear Interpolation")));
432 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
434 values.resize(numInterpolationPoints);
435 for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
436 values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
440 for(LO dof = 0; dof < dofsPerNode; ++dof) {
441 rowIdx = nodeIdx*dofsPerNode + dof;
442 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
443 PCrs->replaceLocalValues(rowIdx, colIndices, values());
449 *out <<
"The calculation of the interpolation stencils has completed." << std::endl;
451 PCrs->fillComplete();
454 *out <<
"All values in P have been set and fillComplete has been performed." << std::endl;
462 if(removeSmallEntries) {
463 *out <<
"Entering remove small entries" << std::endl;
466 ArrayRCP<const size_t> rowptrOrig;
467 ArrayRCP<const LO> colindOrig;
468 ArrayRCP<const Scalar> valuesOrig;
469 PCrs->getAllValues(rowptrOrig, colindOrig, valuesOrig);
471 const size_t numRows =
static_cast<size_t>(rowptrOrig.size() - 1);
472 ArrayRCP<size_t> rowPtr(numRows + 1);
473 ArrayRCP<size_t> nnzOnRows(numRows);
475 size_t countRemovedEntries = 0;
476 for(
size_t rowIdx = 0; rowIdx < numRows; ++rowIdx) {
477 for(
size_t entryIdx = rowptrOrig[rowIdx]; entryIdx < rowptrOrig[rowIdx + 1]; ++entryIdx) {
478 if(Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) < 1e-6) {++countRemovedEntries;}
480 rowPtr[rowIdx + 1] = rowptrOrig[rowIdx + 1] - countRemovedEntries;
481 nnzOnRows[rowIdx] = rowPtr[rowIdx + 1] - rowPtr[rowIdx];
483 GetOStream(
Statistics1) <<
"interp: number of small entries removed= " << countRemovedEntries <<
" / " << rowptrOrig[numRows] << std::endl;
485 size_t countKeptEntries = 0;
486 ArrayRCP<LO> colInd(rowPtr[numRows]);
487 ArrayRCP<SC> values(rowPtr[numRows]);
488 for(
size_t entryIdx = 0; entryIdx < rowptrOrig[numRows]; ++entryIdx) {
489 if(Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) > 1e-6) {
490 colInd[countKeptEntries] = colindOrig[entryIdx];
491 values[countKeptEntries] = valuesOrig[entryIdx];
496 P = rcp(
new CrsMatrixWrap(prolongatorGraph->getRowMap(),
497 prolongatorGraph->getColMap(),
499 RCP<CrsMatrix> PCrsSqueezed = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
500 PCrsSqueezed->resumeFill();
501 PCrsSqueezed->setAllValues(rowPtr, colInd, values);
502 PCrsSqueezed->expertStaticFillComplete(prolongatorGraph->getDomainMap(),
503 prolongatorGraph->getRangeMap());
506 std::vector<size_t> strideInfo(1);
507 strideInfo[0] = dofsPerNode;
508 RCP<const StridedMap> stridedDomainMap =
509 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
511 *out <<
"The strided maps of P have been computed" << std::endl;
514 if (A->IsView(
"stridedMaps") ==
true) {
515 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
517 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
523 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
526 const Array<Array<real_type> > coord,
527 Array<real_type>& stencil)
const {
544 Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
545 Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
546 Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
547 Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
548 Teuchos::SerialDenseSolver<LO,real_type> problem;
549 int iter = 0, max_iter = 5;
550 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
551 paramCoords.size(numDimensions);
553 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
556 solutionDirection.size(numDimensions);
557 residual.size(numDimensions);
561 GetInterpolationFunctions(numDimensions, paramCoords, functions);
562 for(LO i = 0; i < numDimensions; ++i) {
563 residual(i) = coord[0][i];
564 for(LO k = 0; k < numInterpolationPoints; ++k) {
565 residual(i) -= functions[0][k]*coord[k+1][i];
568 norm_ref += residual(i)*residual(i);
569 if(i == numDimensions - 1) {
570 norm_ref = std::sqrt(norm_ref);
574 for(LO j = 0; j < numDimensions; ++j) {
575 for(LO k = 0; k < numInterpolationPoints; ++k) {
576 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
582 problem.setMatrix(Teuchos::rcp(&Jacobian,
false));
583 problem.setVectors(Teuchos::rcp(&solutionDirection,
false), Teuchos::rcp(&residual,
false));
584 if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(
true);}
587 for(LO i = 0; i < numDimensions; ++i) {
588 paramCoords(i) = paramCoords(i) + solutionDirection(i);
592 GetInterpolationFunctions(numDimensions, paramCoords, functions);
593 for(LO i = 0; i < numDimensions; ++i) {
595 for(LO k = 0; k < numInterpolationPoints; ++k) {
596 tmp -= functions[0][k]*coord[k+1][i];
601 norm2 = std::sqrt(norm2);
605 for(LO i = 0; i < numInterpolationPoints; ++i) {
606 stencil[i] = functions[0][i];
611 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
614 const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
616 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
617 if(numDimensions == 1) {
618 xi = parametricCoordinates[0];
620 }
else if(numDimensions == 2) {
621 xi = parametricCoordinates[0];
622 eta = parametricCoordinates[1];
624 }
else if(numDimensions == 3) {
625 xi = parametricCoordinates[0];
626 eta = parametricCoordinates[1];
627 zeta = parametricCoordinates[2];
631 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
632 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
633 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
634 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
635 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
636 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
637 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
638 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
640 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
641 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
642 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
643 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
644 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
645 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
646 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
647 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
649 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
650 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
651 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
652 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
653 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
654 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
655 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
656 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
658 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
659 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
660 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
661 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
662 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
663 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
664 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
665 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
#define SET_VALID_ENTRY(name)
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.
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
void BuildLinearP(Level &coarseLevel, RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, const bool keepD2, RCP< Matrix > &P) const
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
Class that holds all level-specific information.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.