MueLu  Version of the Day
MueLu_GeometricInterpolationPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
48 
49 #include "Xpetra_CrsGraph.hpp"
51 
52 #include "MueLu_MasterList.hpp"
53 #include "MueLu_Monitor.hpp"
54 #include "MueLu_Aggregates.hpp"
55 
56 // Including this one last ensure that the short names of the above headers are defined properly
58 
59 namespace MueLu {
60 
61  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
63  RCP<ParameterList> validParamList = rcp(new ParameterList());
64 
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
66  SET_VALID_ENTRY("gmg: interpolation order");
67  SET_VALID_ENTRY("gmg: build coarse coordinates");
68 #undef SET_VALID_ENTRY
69 
70  // general variables needed in GeometricInterpolationPFactory
71  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
72  "Generating factory of the matrix A");
73  validParamList->set<RCP<const FactoryBase> >("Aggregates", Teuchos::null,
74  "Aggregates generated by StructuredAggregationFactory used to construct a piece-constant prolongator.");
75  validParamList->set<RCP<const FactoryBase> >("prolongatorGraph", Teuchos::null,
76  "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
77  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
78  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
79  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesFineMap", Teuchos::null,
80  "map of the coarse coordinates' GIDs as indexed on the fine mesh.");
81  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesMap", Teuchos::null,
82  "map of the coarse coordinates' GIDs as indexed on the coarse mesh.");
83  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
84  "Fine level nullspace used to construct the coarse level nullspace.");
85  validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
86  "Number of nodes per spatial dimension on the coarse grid.");
87  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
88  "Number of spatial dimensions of the mesh.");
89  validParamList->set<RCP<const FactoryBase> >("interpolationOrder", Teuchos::null,
90  "Interpolation order use to construct the prolongator.");
91 
92  return validParamList;
93  }
94 
95  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
97  DeclareInput(Level& fineLevel, Level& coarseLevel) const {
98  const ParameterList& pL = GetParameterList();
99 
100  Input(fineLevel, "A");
101  Input(fineLevel, "Nullspace");
102  Input(fineLevel, "prolongatorGraph");
103  Input(fineLevel, "lCoarseNodesPerDim");
104  Input(fineLevel, "numDimensions");
105  Input(fineLevel, "interpolationOrder");
106 
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");
112  }
113 
114  }
115 
116  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118  Build(Level& fineLevel, Level &coarseLevel) const {
119  return BuildP(fineLevel, coarseLevel);
120  }
121 
122  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
124  BuildP(Level &fineLevel, Level &coarseLevel) const {
125  FactoryMonitor m(*this, "BuildP", coarseLevel);
126 
127  // Set debug outputs based on environment variable
129  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
130  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
131  out->setShowAllFrontMatter(false).setShowProcRank(true);
132  } else {
133  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
134  }
135 
136  *out << "Starting GeometricInterpolationPFactory::BuildP." << std::endl;
137 
138  // Get inputs from the parameter list
139  const ParameterList& pL = GetParameterList();
140  const bool buildCoarseCoordinates = pL.get<bool>("gmg: build coarse coordinates");
141 
142  // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
143  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
144  RCP<CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
145  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
146  RCP<Matrix> P;
147  const int numDimensions = Get<int>(fineLevel, "numDimensions");
148  const int interpolationOrder = Get<int>(fineLevel, "interpolationOrder");
149 
150  // Check if we need to build coarse coordinates as they are used if we construct
151  // a linear interpolation prolongator
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);
162 
163  Set(coarseLevel, "Coordinates", coarseCoordinates);
164  }
165 
166  *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
167 
168  if(interpolationOrder == 0) {
169  // Compute the prolongator using piece-wise constant interpolation
170  BuildConstantP(P, prolongatorGraph, A);
171  } else if(interpolationOrder == 1) {
172  // Compute the prolongator using piece-wise linear interpolation
173  // First get all the required coordinates to compute the local part of P
174  RCP<realvaluedmultivector_type> ghostCoordinates
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);
180 
181  BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
182  }
183 
184  *out << "The prolongator matrix has been build." << std::endl;
185 
186  // Build the coarse nullspace
187  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
188  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
189  fineNullspace->getNumVectors());
190  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
192  Set(coarseLevel, "Nullspace", coarseNullspace);
193 
194  *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
195 
196  Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
197  Set(coarseLevel, "lNodesPerDim", lNodesPerDir);
198  Set(coarseLevel, "P", P);
199 
200  *out << "GeometricInterpolationPFactory::BuildP has completed." << std::endl;
201 
202  } // BuildP
203 
204  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
206  BuildConstantP(RCP<Matrix>& P, RCP<CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
207 
208  // Set debug outputs based on environment variable
210  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
211  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
212  out->setShowAllFrontMatter(false).setShowProcRank(true);
213  } else {
214  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
215  }
216 
217  // Get the number of dofs per nodes from A and use that number to manually unamalgamate
218  // the maps of the prolongator graph. Eventually this operation will no longer
219  int dofsPerNode = A->GetFixedBlockSize();
220  ArrayView<const GO> initialDomainMapLIDs =
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;
227  }
228  }
229  RCP<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
230  prolongatorGraph->getGlobalNumCols()*dofsPerNode,
231  domainMapLIDs(),
232  prolongatorGraph->getIndexBase(),
233  prolongatorGraph->getComm(),
234  prolongatorGraph->getRowMap()->getNode());
235 
236  ArrayView<const GO> initialColMapLIDs =
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;
243  }
244  }
245  RCP<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
246  prolongatorGraph->getGlobalNumCols()*dofsPerNode,
247  colMapLIDs(),
248  prolongatorGraph->getIndexBase(),
249  prolongatorGraph->getComm(),
250  prolongatorGraph->getColMap()->getNode());
251 
252  std::vector<size_t> strideInfo(1);
253  strideInfo[0] = dofsPerNode;
254  RCP<const StridedMap> stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo);
255 
256  // Create the prolongator matrix and its associated objects
257  P = rcp(new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile));
258  // P = rcp(new CrsMatrixWrap(graph));
259  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
260 
261  ArrayRCP<size_t> iaP;
262  ArrayRCP<LO> jaP;
263  ArrayRCP<SC> valP;
264 
265  PCrs->allocateAllValues(A->getDomainMap()->getNodeNumElements(), iaP, jaP, valP);
266 
267  ArrayView<size_t> ia = iaP();
268  ArrayView<LO> ja = jaP();
269  ArrayView<SC> val = valP();
270  ia[0] = 0;
271 
272  // Actually fill up the matrix, in this case all values are one, pretty easy!
273  int dofIdx;
274  ArrayView<const LO> colIdx;
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;
281  val[dofIdx] = 1.0;
282  }
283  }
284 
285  // call fill complete on the prolongator
286  PCrs->setAllValues(iaP, jaP, valP);
287  PCrs->expertStaticFillComplete(domainMap, A->getDomainMap());
288 
289  // set StridingInformation of P
290  if (A->IsView("stridedMaps") == true) {
291  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
292  } else {
293  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
294  }
295 
296  } // BuildConstantP
297 
298  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
301  RCP<realvaluedmultivector_type>& fineCoordinates,
302  RCP<realvaluedmultivector_type>& ghostCoordinates,
303  const int numDimensions, RCP<Matrix>& P) const {
304 
305  // Set debug outputs based on environment variable
307  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
308  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
309  out->setShowAllFrontMatter(false).setShowProcRank(true);
310  } else {
311  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
312  }
313 
314  *out << "Entering BuildLinearP" << std::endl;
315 
316  // Extract coordinates for interpolation stencil calculations
317  Array<ArrayView<const real_type> > fineCoords(3);
318  Array<ArrayView<const real_type> > ghostCoords(3);
319  real_type realZero = Teuchos::as<real_type>(0.0);
320  Array<real_type> fineZero(fineCoordinates->getLocalLength(), realZero);
321  Array<real_type> ghostZero(ghostCoordinates->getLocalLength(), realZero);
322  for(int dim = 0; dim < 3; ++dim) {
323  if(dim < numDimensions) {
324  fineCoords[dim] = fineCoordinates->getData(dim)();
325  ghostCoords[dim] = ghostCoordinates->getData(dim)();
326  } else {
327  fineCoords[dim] = fineZero();
328  ghostCoords[dim] = ghostZero();
329  }
330  }
331 
332  *out << "Coordinates extracted from the multivectors!" << std::endl;
333 
334  // Compute 2^numDimensions using bit logic to avoid round-off errors
335  const int numInterpolationPoints = 1 << numDimensions;
336  const int dofsPerNode = A->GetFixedBlockSize();
337  ArrayView<const GO> initialDomainMapLIDs =
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;
344  }
345  }
346  RCP<Map> domainMap = MapFactory::Build(prolongatorGraph->getRowMap()->lib(),
347  prolongatorGraph->getGlobalNumCols()*dofsPerNode,
348  domainMapLIDs(),
349  prolongatorGraph->getIndexBase(),
350  prolongatorGraph->getComm(),
351  prolongatorGraph->getRowMap()->getNode());
352 
353  ArrayView<const GO> initialColMapLIDs =
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;
360  }
361  }
362  RCP<Map> colMap = MapFactory::Build(prolongatorGraph->getColMap()->lib(),
363  prolongatorGraph->getGlobalNumCols()*dofsPerNode,
364  colMapLIDs(),
365  prolongatorGraph->getIndexBase(),
366  prolongatorGraph->getComm(),
367  prolongatorGraph->getColMap()->getNode());
368 
369  std::vector<size_t> strideInfo(1);
370  strideInfo[0] = dofsPerNode;
371  RCP<const StridedMap> stridedDomainMap = StridedMapFactory::Build(domainMap, strideInfo);
372 
373  *out << "The maps of P have been computed" << std::endl;
374 
375  P = rcp(new CrsMatrixWrap(A->getDomainMap(), colMap, 0, Xpetra::StaticProfile));
376  // P = rcp(new CrsMatrixWrap(graph));
377  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
378 
379  ArrayRCP<size_t> iaP;
380  ArrayRCP<LO> jaP;
381  ArrayRCP<SC> valP;
382 
383  LO numNonZeroP = prolongatorGraph->getNodeNumEntries()*dofsPerNode;
384  PCrs->allocateAllValues(numNonZeroP, iaP, jaP, valP);
385 
386  *out << "dofsPerNode=" << dofsPerNode << std::endl;
387  *out << "number of non-zeroes in P: " << numNonZeroP << std::endl;
388 
389  ArrayView<size_t> ia = iaP();
390  ArrayView<LO> ja = jaP();
391  ArrayView<SC> val = valP();
392  ia[0] = 0;
393 
394  LO rowDofIdx, colDofIdx;
395  ArrayView<const LO> colIndices;
396  Array<Array<real_type> > coords(numInterpolationPoints + 1);
397  Array<real_type> stencil(numInterpolationPoints);
398  for(LO rowIdx = 0; rowIdx < prolongatorGraph->getNodeNumRows(); ++rowIdx) {
399  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
400 
401  // rowIdx and colIdx correspond to single ids of nodes on the fine resp. coarse mesh
402  // rowDofIdx and colDofIdx correspond to single ids of dofs on the fine resp. coarse mesh
403  if(colIndices.size() == 1) {
404  // If colIndices.size() == 1, we are handling a coarse node
405  // we only need to stick a one in the correct location!
406  for(int dof = 0; dof < dofsPerNode; ++dof) {
407  rowDofIdx = rowIdx*dofsPerNode + dof;
408  ia[rowDofIdx + 1] = ia[rowDofIdx] + 1;
409 
410  colDofIdx = ia[rowDofIdx];
411  ja[colDofIdx] = colIndices[0]*dofsPerNode + dof;
412  val[colDofIdx] = 1.0;
413  }
414  } else {
415  // If colIndices.size() > 1, we need to compute the linear interpolation coefficients
416  for(int dof = 0; dof < dofsPerNode; ++dof) {
417  rowDofIdx = rowIdx*dofsPerNode + dof;
418  ia[rowDofIdx + 1] = ia[rowDofIdx] + colIndices.size();
419 
420  // Extract the coordinates associated with the current node
421  // and the neighboring coarse nodes
422  coords[0].resize(3);
423  for(int dim = 0; dim < 3; ++dim) {
424  coords[0][dim] = fineCoords[dim][rowIdx];
425  }
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]];
430  }
431  }
432  ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
433 
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];
438  }
439  }
440  }
441  }
442 
443  *out << "The calculation of the interpolation stencils has completed." << std::endl;
444 
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());
448 
449  *out << "All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
450 
451  // set StridingInformation of P
452  if (A->IsView("stridedMaps") == true) {
453  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
454  } else {
455  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
456  }
457 
458  } // BuildLinearP
459 
460 
461  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
463  ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
464  const Array<Array<real_type> > coord,
465  Array<real_type>& stencil) const {
466 
467  // 7 8 Find xi, eta and zeta such that
468  // x---------x
469  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
470  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
471  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
472  // | | *P | |
473  // | x------|--x We can do this with a Newton solver:
474  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
475  // |/ |/ We compute the Jacobian and iterate until convergence...
476  // z y x---------x
477  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
478  // |/ give us the weights for the interpolation stencil!
479  // o---x
480  //
481 
482  Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
483  Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
484  Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
485  Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
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);
490 
491  while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
492  ++iter;
493  norm2 = 0.0;
494  solutionDirection.size(numDimensions);
495  residual.size(numDimensions);
496  Jacobian = 0.0;
497 
498  // Compute Jacobian and Residual
499  GetInterpolationFunctions(numDimensions, paramCoords, functions);
500  for(LO i = 0; i < numDimensions; ++i) {
501  residual(i) = coord[0][i]; // Add coordinates from point of interest
502  for(LO k = 0; k < numInterpolationPoints; ++k) {
503  residual(i) -= functions[0][k]*coord[k+1][i]; //Remove contribution from all coarse points
504  }
505  if(iter == 1) {
506  norm_ref += residual(i)*residual(i);
507  if(i == numDimensions - 1) {
508  norm_ref = std::sqrt(norm_ref);
509  }
510  }
511 
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];
515  }
516  }
517  }
518 
519  // Set Jacobian, Vectors and solve problem
520  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
521  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
522  if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(true);}
523  problem.solve();
524 
525  for(LO i = 0; i < numDimensions; ++i) {
526  paramCoords(i) = paramCoords(i) + solutionDirection(i);
527  }
528 
529  // Recompute Residual norm
530  GetInterpolationFunctions(numDimensions, paramCoords, functions);
531  for(LO i = 0; i < numDimensions; ++i) {
532  real_type tmp = coord[0][i];
533  for(LO k = 0; k < numInterpolationPoints; ++k) {
534  tmp -= functions[0][k]*coord[k+1][i];
535  }
536  norm2 += tmp*tmp;
537  tmp = 0.0;
538  }
539  norm2 = std::sqrt(norm2);
540  }
541 
542  // Load the interpolation values onto the stencil.
543  for(LO i = 0; i < 8; ++i) {
544  stencil[i] = functions[0][i];
545  }
546 
547  } // End ComputeLinearInterpolationStencil
548 
549  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
551  GetInterpolationFunctions(const LO numDimensions,
552  const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
553  real_type functions[4][8]) const {
554  real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
555  if(numDimensions == 1) {
556  xi = parametricCoordinates[0];
557  denominator = 2.0;
558  } else if(numDimensions == 2) {
559  xi = parametricCoordinates[0];
560  eta = parametricCoordinates[1];
561  denominator = 4.0;
562  } else if(numDimensions == 3) {
563  xi = parametricCoordinates[0];
564  eta = parametricCoordinates[1];
565  zeta = parametricCoordinates[2];
566  denominator = 8.0;
567  }
568 
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;
577 
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;
586 
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;
595 
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;
604 
605  } // End GetInterpolationFunctions
606 
607 } // namespace MueLu
608 
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.
LocalOrdinal LO
size_type size() const
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.
Definition: MueLu_Level.hpp:99
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
static void sortCrsEntries(const Teuchos::ArrayView< size_t > &CRS_rowptr, const Teuchos::ArrayView< LocalOrdinal > &CRS_colind, const Teuchos::ArrayView< Scalar > &CRS_vals, const UnderlyingLib lib)
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
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
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