46 #ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_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 spacial dimensions in the problem.");
88 "Number of nodes per spatial dimension on the coarse grid.");
90 return validParamList;
93 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 Input(fineLevel,
"A");
99 Input(fineLevel,
"Nullspace");
100 Input(fineLevel,
"numDimensions");
101 Input(fineLevel,
"prolongatorGraph");
102 Input(fineLevel,
"lCoarseNodesPerDim");
104 if( pL.
get<
bool>(
"interp: build coarse coordinates") ||
105 (pL.
get<
int>(
"interp: interpolation order") == 1) ) {
106 Input(fineLevel,
"Coordinates");
107 Input(fineLevel,
"coarseCoordinatesFineMap");
108 Input(fineLevel,
"coarseCoordinatesMap");
113 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 return BuildP(fineLevel, coarseLevel);
119 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
127 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
133 *out <<
"Starting GeometricInterpolationPFactory::BuildP." << std::endl;
137 const bool buildCoarseCoordinates = pL.
get<
bool>(
"interp: build coarse coordinates");
138 const int interpolationOrder = pL.
get<
int> (
"interp: interpolation order");
139 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
149 if(buildCoarseCoordinates || (interpolationOrder == 1)) {
151 RCP<const Map> coarseCoordsFineMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesFineMap");
152 RCP<const Map> coarseCoordsMap = Get< RCP<const Map> >(fineLevel,
"coarseCoordinatesMap");
153 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
155 fineCoordinates->getNumVectors());
156 RCP<const Import> coordsImporter = ImportFactory::Build(fineCoordinates->getMap(),
157 coarseCoordsFineMap);
158 coarseCoordinates->doImport(*fineCoordinates, *coordsImporter,
Xpetra::INSERT);
159 coarseCoordinates->replaceMap(coarseCoordsMap);
161 Set(coarseLevel,
"Coordinates", coarseCoordinates);
164 *out <<
"Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
166 if(interpolationOrder == 0) {
169 BuildConstantP(P, prolongatorGraph, A);
170 }
else if(interpolationOrder == 1) {
175 fineCoordinates->getNumVectors());
176 RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
177 prolongatorGraph->getColMap());
178 ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter,
Xpetra::INSERT);
181 BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
184 *out <<
"The prolongator matrix has been built." << std::endl;
189 RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel,
"Nullspace");
190 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
191 fineNullspace->getNumVectors());
194 Set(coarseLevel,
"Nullspace", coarseNullspace);
197 *out <<
"The coarse nullspace is constructed and set on the coarse level." << std::endl;
199 Set(coarseLevel,
"P", P);
201 *out <<
"GeometricInterpolationPFactory::BuildP has completed." << std::endl;
205 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
211 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
212 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
218 *out <<
"BuildConstantP" << std::endl;
220 std::vector<size_t> strideInfo(1);
221 strideInfo[0] = A->GetFixedBlockSize();
223 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
225 *out <<
"Call prolongator constructor" << std::endl;
229 P =
rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
230 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
231 PCrs->setAllToScalar(1.0);
232 PCrs->fillComplete();
235 if (A->IsView(
"stridedMaps") ==
true) {
236 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
238 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
243 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
252 if(
const char* dbg = std::getenv(
"MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
253 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
259 *out <<
"Entering BuildLinearP" << std::endl;
262 const LO numFineNodes = fineCoordinates->getLocalLength();
263 const LO numGhostNodes = ghostCoordinates->getLocalLength();
266 const real_type realZero = Teuchos::as<real_type>(0.0);
269 for(
int dim = 0; dim < 3; ++dim) {
270 if(dim < numDimensions) {
271 fineCoords[dim] = fineCoordinates->getData(dim);
272 ghostCoords[dim] = ghostCoordinates->getData(dim);
274 fineCoords[dim] = fineZero;
275 ghostCoords[dim] = ghostZero;
279 *out <<
"Coordinates extracted from the multivectors!" << std::endl;
282 const int numInterpolationPoints = 1 << numDimensions;
283 const int dofsPerNode = A->GetFixedBlockSize();
285 std::vector<size_t> strideInfo(1);
286 strideInfo[0] = dofsPerNode;
288 StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
290 *out <<
"The maps of P have been computed" << std::endl;
293 P =
rcp(
new CrsMatrixWrap(prolongatorGraph, dummyList));
294 RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
297 LO interpolationNodeIdx = 0, rowIdx = 0;
302 for(
LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
303 if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
306 for(
LO dof = 0; dof < dofsPerNode; ++dof) {
307 rowIdx = nodeIdx*dofsPerNode + dof;
308 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
309 PCrs->replaceLocalValues(rowIdx, colIndices, values());
315 for(
int dim = 0; dim < 3; ++dim) {
316 coords[0][dim] = fineCoords[dim][nodeIdx];
318 prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
319 for(
int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
320 coords[interpolationIdx + 1].
resize(3);
321 interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
322 for(
int dim = 0; dim < 3; ++dim) {
323 coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
326 ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
327 values.
resize(numInterpolationPoints);
328 for(
LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
329 values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
333 for(
LO dof = 0; dof < dofsPerNode; ++dof) {
334 rowIdx = nodeIdx*dofsPerNode + dof;
335 prolongatorGraph->getLocalRowView(rowIdx, colIndices);
336 PCrs->replaceLocalValues(rowIdx, colIndices, values());
341 *out <<
"The calculation of the interpolation stencils has completed." << std::endl;
343 PCrs->fillComplete();
345 *out <<
"All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
348 if (A->IsView(
"stridedMaps") ==
true) {
349 P->CreateView(
"stridedMaps", A->getRowMap(
"stridedMaps"), stridedDomainMap);
351 P->CreateView(
"stridedMaps", P->getRangeMap(), stridedDomainMap);
357 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
383 int iter = 0, max_iter = 5;
384 real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
385 paramCoords.
size(numDimensions);
387 while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
390 solutionDirection.
size(numDimensions);
391 residual.
size(numDimensions);
395 GetInterpolationFunctions(numDimensions, paramCoords, functions);
396 for(
LO i = 0; i < numDimensions; ++i) {
397 residual(i) = coord[0][i];
398 for(
LO k = 0; k < numInterpolationPoints; ++k) {
399 residual(i) -= functions[0][k]*coord[k+1][i];
402 norm_ref += residual(i)*residual(i);
403 if(i == numDimensions - 1) {
404 norm_ref = std::sqrt(norm_ref);
408 for(
LO j = 0; j < numDimensions; ++j) {
409 for(
LO k = 0; k < numInterpolationPoints; ++k) {
410 Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
421 for(
LO i = 0; i < numDimensions; ++i) {
422 paramCoords(i) = paramCoords(i) + solutionDirection(i);
426 GetInterpolationFunctions(numDimensions, paramCoords, functions);
427 for(
LO i = 0; i < numDimensions; ++i) {
429 for(
LO k = 0; k < numInterpolationPoints; ++k) {
430 tmp -= functions[0][k]*coord[k+1][i];
435 norm2 = std::sqrt(norm2);
439 for(
LO i = 0; i < numInterpolationPoints; ++i) {
440 stencil[i] = functions[0][i];
445 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
450 real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
451 if(numDimensions == 1) {
452 xi = parametricCoordinates[0];
454 }
else if(numDimensions == 2) {
455 xi = parametricCoordinates[0];
456 eta = parametricCoordinates[1];
458 }
else if(numDimensions == 3) {
459 xi = parametricCoordinates[0];
460 eta = parametricCoordinates[1];
461 zeta = parametricCoordinates[2];
465 functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
466 functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
467 functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
468 functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
469 functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
470 functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
471 functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
472 functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
474 functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
475 functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
476 functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
477 functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
478 functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
479 functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
480 functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
481 functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
483 functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
484 functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
485 functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
486 functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
487 functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
488 functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
489 functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
490 functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
492 functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
493 functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
494 functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
495 functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
496 functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
497 functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
498 functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
499 functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
505 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
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)
Timer to be used in factories. Similar to Monitor but with additional timers.
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Class that holds all level-specific information.
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void factorWithEquilibration(bool flag)
#define SET_VALID_ENTRY(name)
void BuildLinearP(RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, RCP< Matrix > &P) const
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
void resize(size_type new_size, const value_type &x=value_type())
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
int size(OrdinalType length_in)
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)