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" 61 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 68 #undef SET_VALID_ENTRY 72 "Generating factory of the matrix A");
74 "Aggregates generated by StructuredAggregationFactory used to construct a piece-constant prolongator.");
76 "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
78 "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
80 "map of the coarse coordinates' GIDs as indexed on the fine mesh.");
82 "map of the coarse coordinates' GIDs as indexed on the coarse mesh.");
84 "Fine level nullspace used to construct the coarse level nullspace.");
86 "Number of nodes per spatial dimension on the coarse grid.");
88 "Number of spatial dimensions of the mesh.");
90 "Interpolation order use to construct the prolongator.");
92 return validParamList;
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 Input(fineLevel,
"A");
101 Input(fineLevel,
"Nullspace");
102 Input(fineLevel,
"prolongatorGraph");
103 Input(fineLevel,
"lCoarseNodesPerDim");
104 Input(fineLevel,
"numDimensions");
105 Input(fineLevel,
"interpolationOrder");
107 if( pL.
get<
bool>(
"gmg: build coarse coordinates") ||
108 (pL.
get<
int>(
"gmg: interpolation order") == 1) ) {
109 Input(fineLevel,
"Coordinates");
110 Input(fineLevel,
"coarseCoordinatesFineMap");
111 Input(fineLevel,
"coarseCoordinatesMap");
116 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
119 return BuildP(fineLevel, coarseLevel);
122 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
129 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
130 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
136 *out <<
"Starting GeometricInterpolationPFactory::BuildP." << std::endl;
140 const bool buildCoarseCoordinates = pL.
get<
bool>(
"gmg: build coarse coordinates");
144 RCP<CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel,
"prolongatorGraph");
147 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
148 const int interpolationOrder = Get<int>(fineLevel,
"interpolationOrder");
152 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
153 RCP<const Map> coarseCoordsFineMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesFineMap");
154 RCP<const Map> coarseCoordsMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesMap");
155 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
156 coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::Build(coarseCoordsFineMap,
157 fineCoordinates->getNumVectors());
158 RCP<const Import> coordsImporter = ImportFactory::Build(fineCoordinates->getMap(),
159 coarseCoordsFineMap);
160 coarseCoordinates->doImport(*fineCoordinates, *coordsImporter, Xpetra::INSERT);
161 coarseCoordinates->replaceMap(coarseCoordsMap);
163 Set(coarseLevel,
"Coordinates", coarseCoordinates);
166 *out <<
"Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
168 if(interpolationOrder == 0) {
170 BuildConstantP(P, prolongatorGraph, A);
171 }
else if(interpolationOrder == 1) {
175 = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(prolongatorGraph->getColMap(),
176 fineCoordinates->getNumVectors());
177 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
178 prolongatorGraph->getColMap());
179 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
181 BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
184 *out <<
"The prolongator matrix has been build." << std::endl;
187 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
188 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
189 fineNullspace->getNumVectors());
192 Set(coarseLevel,
"Nullspace", coarseNullspace);
194 *out <<
"The coarse nullspace is constructed and set on the coarse level." << std::endl;
196 Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel,
"lCoarseNodesPerDim");
197 Set(coarseLevel,
"lNodesPerDim", lNodesPerDir);
198 Set(coarseLevel,
"P", P);
200 *out <<
"GeometricInterpolationPFactory::BuildP has completed." << std::endl;
204 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
210 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
211 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
219 int dofsPerNode = A->GetFixedBlockSize();
221 prolongatorGraph->getDomainMap()->getNodeElementList();
222 Array<GO> domainMapLIDs(initialDomainMapLIDs.
size()*dofsPerNode);
223 for(LO elementIdx = 0; elementIdx < as<LO>(domainMapLIDs.size()); ++elementIdx) {
224 for(
int dof = 0; dof < dofsPerNode; ++dof) {
225 domainMapLIDs[elementIdx*dofsPerNode + dof] =
226 initialDomainMapLIDs[elementIdx]*dofsPerNode + dof;
229 RCP<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
230 prolongatorGraph->getGlobalNumCols()*dofsPerNode,
232 prolongatorGraph->getIndexBase(),
233 prolongatorGraph->getComm(),
234 prolongatorGraph->getRowMap()->getNode());
237 prolongatorGraph->getColMap()->getNodeElementList();
238 Array<GO> colMapLIDs(initialColMapLIDs.size()*dofsPerNode);
239 for(LO elementIdx = 0; elementIdx < as<LO>(colMapLIDs.size()); ++elementIdx) {
240 for(
int dof = 0; dof < dofsPerNode; ++dof) {
241 colMapLIDs[elementIdx*dofsPerNode + dof] =
242 initialColMapLIDs[elementIdx]*dofsPerNode + dof;
245 RCP<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
246 prolongatorGraph->getGlobalNumCols()*dofsPerNode,
248 prolongatorGraph->getIndexBase(),
249 prolongatorGraph->getComm(),
250 prolongatorGraph->getColMap()->getNode());
252 std::vector<size_t> strideInfo(1);
253 strideInfo[0] = dofsPerNode;
257 P =
rcp(
new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile));
259 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
265 PCrs->allocateAllValues(A->getDomainMap()->getNodeNumElements(), iaP, jaP, valP);
275 for(LO rowIdx = 0; rowIdx < prolongatorGraph->getNodeNumRows(); ++rowIdx) {
276 prolongatorGraph->getLocalRowView(rowIdx, colIdx);
277 for(
int dof = 0; dof < dofsPerNode; ++dof) {
278 dofIdx = rowIdx*dofsPerNode + dof;
279 ia[dofIdx + 1] = dofIdx + 1;
280 ja[dofIdx] = colIdx[0]*dofsPerNode + dof;
286 PCrs->setAllValues(iaP, jaP, valP);
287 PCrs->expertStaticFillComplete(domainMap, A->getDomainMap());
290 if (A->IsView(
"stridedMaps") ==
true) {
291 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
293 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
298 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
307 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
308 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
314 *out <<
"Entering BuildLinearP" << std::endl;
319 real_type realZero = Teuchos::as<real_type>(0.0);
322 for(
int dim = 0; dim < 3; ++dim) {
323 if(dim < numDimensions) {
324 fineCoords[dim] = fineCoordinates->getData(dim)();
325 ghostCoords[dim] = ghostCoordinates->getData(dim)();
327 fineCoords[dim] = fineZero();
328 ghostCoords[dim] = ghostZero();
332 *out <<
"Coordinates extracted from the multivectors!" << std::endl;
335 const int numInterpolationPoints = 1 << numDimensions;
336 const int dofsPerNode = A->GetFixedBlockSize();
338 prolongatorGraph->getDomainMap()->getNodeElementList();
339 Array<GO> domainMapLIDs(initialDomainMapLIDs.
size()*dofsPerNode);
340 for(LO elementIdx = 0; elementIdx < as<LO>(domainMapLIDs.size()); ++elementIdx) {
341 for(
int dof = 0; dof < dofsPerNode; ++dof) {
342 domainMapLIDs[elementIdx*dofsPerNode + dof] =
343 initialDomainMapLIDs[elementIdx]*dofsPerNode + dof;
346 RCP<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
347 prolongatorGraph->getGlobalNumCols()*dofsPerNode,
349 prolongatorGraph->getIndexBase(),
350 prolongatorGraph->getComm(),
351 prolongatorGraph->getRowMap()->getNode());
354 prolongatorGraph->getColMap()->getNodeElementList();
355 Array<GO> colMapLIDs(initialColMapLIDs.size()*dofsPerNode);
356 for(LO elementIdx = 0; elementIdx < as<LO>(colMapLIDs.size()); ++elementIdx) {
357 for(
int dof = 0; dof < dofsPerNode; ++dof) {
358 colMapLIDs[elementIdx*dofsPerNode + dof] =
359 initialColMapLIDs[elementIdx]*dofsPerNode + dof;
362 RCP<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
363 prolongatorGraph->getGlobalNumCols()*dofsPerNode,
365 prolongatorGraph->getIndexBase(),
366 prolongatorGraph->getComm(),
367 prolongatorGraph->getColMap()->getNode());
369 std::vector<size_t> strideInfo(1);
370 strideInfo[0] = dofsPerNode;
373 *out <<
"The maps of P have been computed" << std::endl;
375 P =
rcp(
new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile));
377 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
383 LO numNonZeroP = prolongatorGraph->getNodeNumEntries()*dofsPerNode;
384 PCrs->allocateAllValues(numNonZeroP, iaP, jaP, valP);
386 *out <<
"dofsPerNode=" << dofsPerNode << std::endl;
387 *out <<
"number of non-zeroes in P: " << numNonZeroP << std::endl;
394 LO rowDofIdx, colDofIdx;
398 for(LO rowIdx = 0; rowIdx < prolongatorGraph->getNodeNumRows(); ++rowIdx) {
399 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
403 if(colIndices.
size() == 1) {
406 for(
int dof = 0; dof < dofsPerNode; ++dof) {
407 rowDofIdx = rowIdx*dofsPerNode + dof;
408 ia[rowDofIdx + 1] = ia[rowDofIdx] + 1;
410 colDofIdx = ia[rowDofIdx];
411 ja[colDofIdx] = colIndices[0]*dofsPerNode + dof;
412 val[colDofIdx] = 1.0;
416 for(
int dof = 0; dof < dofsPerNode; ++dof) {
417 rowDofIdx = rowIdx*dofsPerNode + dof;
418 ia[rowDofIdx + 1] = ia[rowDofIdx] + colIndices.
size();
423 for(
int dim = 0; dim < 3; ++dim) {
424 coords[0][dim] = fineCoords[dim][rowIdx];
426 for(
int interpolationIdx = 0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
427 coords[interpolationIdx + 1].
resize(3);
428 for(
int dim = 0; dim < 3; ++dim) {
429 coords[interpolationIdx + 1][dim] = ghostCoords[dim][colIndices[interpolationIdx]];
432 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
434 for(
int colIdx = 0; colIdx < colIndices.
size(); ++colIdx) {
435 colDofIdx = ia[rowDofIdx] + colIdx;
436 ja[colDofIdx] = colIndices[colIdx]*dofsPerNode + dof;
437 val[colDofIdx] = stencil[colIdx];
443 *out <<
"The calculation of the interpolation stencils has completed." << std::endl;
445 Xpetra::CrsMatrixUtils<SC,LO,GO,NO>::sortCrsEntries(ia, ja, val, A->getDomainMap()->lib());
446 PCrs->setAllValues(iaP, jaP, valP);
447 PCrs->expertStaticFillComplete(domainMap, A->getDomainMap());
449 *out <<
"All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
452 if (A->IsView(
"stridedMaps") ==
true) {
453 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
455 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
461 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
487 int iter = 0, max_iter = 5;
488 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
489 paramCoords.
size(numDimensions);
491 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
494 solutionDirection.
size(numDimensions);
495 residual.
size(numDimensions);
499 GetInterpolationFunctions(numDimensions, paramCoords, functions);
500 for(LO i = 0; i < numDimensions; ++i) {
501 residual(i) = coord[0][i];
502 for(LO k = 0; k < numInterpolationPoints; ++k) {
503 residual(i) -= functions[0][k]*coord[k+1][i];
506 norm_ref += residual(i)*residual(i);
507 if(i == numDimensions - 1) {
508 norm_ref = std::sqrt(norm_ref);
512 for(LO j = 0; j < numDimensions; ++j) {
513 for(LO k = 0; k < numInterpolationPoints; ++k) {
514 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
525 for(LO i = 0; i < numDimensions; ++i) {
526 paramCoords(i) = paramCoords(i) + solutionDirection(i);
530 GetInterpolationFunctions(numDimensions, paramCoords, functions);
531 for(LO i = 0; i < numDimensions; ++i) {
533 for(LO k = 0; k < numInterpolationPoints; ++k) {
534 tmp -= functions[0][k]*coord[k+1][i];
539 norm2 = std::sqrt(norm2);
543 for(LO i = 0; i < 8; ++i) {
544 stencil[i] = functions[0][i];
549 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
554 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
555 if(numDimensions == 1) {
556 xi = parametricCoordinates[0];
558 }
else if(numDimensions == 2) {
559 xi = parametricCoordinates[0];
560 eta = parametricCoordinates[1];
562 }
else if(numDimensions == 3) {
563 xi = parametricCoordinates[0];
564 eta = parametricCoordinates[1];
565 zeta = parametricCoordinates[2];
569 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
570 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
571 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
572 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
573 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
574 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
575 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
576 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
578 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
579 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
580 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
581 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
582 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
583 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
584 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
585 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
587 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
588 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
589 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
590 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
591 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
592 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
593 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
594 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
596 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
597 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
598 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
599 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
600 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
601 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
602 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
603 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
609 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
typename Teuchos::ScalarTraits< SC >::magnitudeType real_type
basic_FancyOStream & setShowProcRank(const bool showProcRank)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
void factorWithEquilibration(bool flag)
#define SET_VALID_ENTRY(name)
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
void resize(size_type new_size, const value_type &x=value_type())
void BuildLinearP(RCP< Matrix > &A, RCP< CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, RCP< Matrix > &P) const
int size(OrdinalType length_in)
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
void BuildConstantP(RCP< Matrix > &P, RCP< CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const