MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
50 #include "Xpetra_CrsMatrixUtils.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
57 #include "MueLu_GeometricInterpolationPFactory_decl.hpp"
58 
59 namespace MueLu {
60 
61  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
62  RCP<const ParameterList> GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
63  RCP<ParameterList> validParamList = rcp(new ParameterList());
64 
65 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
66  SET_VALID_ENTRY("interp: interpolation order");
67  SET_VALID_ENTRY("interp: 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> >("numDimensions", Teuchos::null,
86  "Number of spacial dimensions in the problem.");
87  validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
88  "Number of nodes per spatial dimension on the coarse grid.");
89 
90  return validParamList;
91  }
92 
93  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
94  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
95  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
96  const ParameterList& pL = GetParameterList();
97 
98  Input(fineLevel, "A");
99  Input(fineLevel, "Nullspace");
100  Input(fineLevel, "numDimensions");
101  Input(fineLevel, "prolongatorGraph");
102  Input(fineLevel, "lCoarseNodesPerDim");
103 
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");
109  }
110 
111  }
112 
113  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
114  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
115  Build(Level& fineLevel, Level &coarseLevel) const {
116  return BuildP(fineLevel, coarseLevel);
117  }
118 
119  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
120  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
121  BuildP(Level &fineLevel, Level &coarseLevel) const {
122  FactoryMonitor m(*this, "BuildP", coarseLevel);
123 
124  // Set debug outputs based on environment variable
125  RCP<Teuchos::FancyOStream> out;
126  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
127  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
128  out->setShowAllFrontMatter(false).setShowProcRank(true);
129  } else {
130  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
131  }
132 
133  *out << "Starting GeometricInterpolationPFactory::BuildP." << std::endl;
134 
135  // Get inputs from the parameter list
136  const ParameterList& pL = GetParameterList();
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");
140 
141  // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
142  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
143  RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
144  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
145  RCP<Matrix> P;
146 
147  // Check if we need to build coarse coordinates as they are used if we construct
148  // a linear interpolation prolongator
149  if(buildCoarseCoordinates || (interpolationOrder == 1)) {
150  SubFactoryMonitor sfm(*this, "BuildCoordinates", coarseLevel);
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");
154  coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::Build(coarseCoordsFineMap,
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);
160 
161  Set(coarseLevel, "Coordinates", coarseCoordinates);
162  }
163 
164  *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
165 
166  if(interpolationOrder == 0) {
167  SubFactoryMonitor sfm(*this, "BuildConstantP", coarseLevel);
168  // Compute the prolongator using piece-wise constant interpolation
169  BuildConstantP(P, prolongatorGraph, A);
170  } else if(interpolationOrder == 1) {
171  // Compute the prolongator using piece-wise linear interpolation
172  // First get all the required coordinates to compute the local part of P
173  RCP<const Map> prolongatorColMap = prolongatorGraph->getColMap();
174 
175  const size_t dofsPerNode = static_cast<size_t>(A->GetFixedBlockSize());
176  const size_t numColIndices = prolongatorColMap->getNodeNumElements();
177  TEUCHOS_TEST_FOR_EXCEPTION((numColIndices % dofsPerNode) != 0,
178  Exceptions::RuntimeError,
179  "Something went wrong, the number of columns in the prolongator is not a multiple of dofsPerNode!");
180  const size_t numGhostCoords = numColIndices / dofsPerNode;
181  const GO indexBase = prolongatorColMap->getIndexBase();
182  const GO coordIndexBase = fineCoordinates->getMap()->getIndexBase();
183 
184  ArrayView<const GO> prolongatorColIndices = prolongatorColMap->getNodeElementList();
185  Array<GO> ghostCoordIndices(numGhostCoords);
186  for(size_t ghostCoordIdx = 0; ghostCoordIdx < numGhostCoords; ++ghostCoordIdx) {
187  ghostCoordIndices[ghostCoordIdx]
188  = (prolongatorColIndices[ghostCoordIdx*dofsPerNode] - indexBase) / dofsPerNode
189  + coordIndexBase;
190  }
191  RCP<Map> ghostCoordMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
192  prolongatorColMap->getGlobalNumElements() / dofsPerNode,
193  ghostCoordIndices(),
194  coordIndexBase,
195  fineCoordinates->getMap()->getComm());
196 
197  RCP<realvaluedmultivector_type> ghostCoordinates
198  = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(ghostCoordMap,
199  fineCoordinates->getNumVectors());
200  RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
201  ghostCoordMap);
202  ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
203 
204  SubFactoryMonitor sfm(*this, "BuildLinearP", coarseLevel);
205  BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
206  }
207 
208  *out << "The prolongator matrix has been built." << std::endl;
209 
210  {
211  SubFactoryMonitor sfm(*this, "BuildNullspace", coarseLevel);
212  // Build the coarse nullspace
213  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
214  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
215  fineNullspace->getNumVectors());
216  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
217  Teuchos::ScalarTraits<SC>::zero());
218  Set(coarseLevel, "Nullspace", coarseNullspace);
219  }
220 
221  *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
222 
223  Set(coarseLevel, "P", P);
224 
225  *out << "GeometricInterpolationPFactory::BuildP has completed." << std::endl;
226 
227  } // BuildP
228 
229  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
230  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
231  BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
232 
233  // Set debug outputs based on environment variable
234  RCP<Teuchos::FancyOStream> out;
235  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
236  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
237  out->setShowAllFrontMatter(false).setShowProcRank(true);
238  } else {
239  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
240  }
241 
242  *out << "BuildConstantP" << std::endl;
243 
244  std::vector<size_t> strideInfo(1);
245  strideInfo[0] = A->GetFixedBlockSize();
246  RCP<const StridedMap> stridedDomainMap =
247  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
248 
249  *out << "Call prolongator constructor" << std::endl;
250 
251  // Create the prolongator matrix and its associated objects
252  RCP<ParameterList> dummyList = rcp(new ParameterList());
253  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
254  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
255  PCrs->setAllToScalar(1.0);
256  PCrs->fillComplete();
257 
258  // set StridingInformation of P
259  if (A->IsView("stridedMaps") == true) {
260  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
261  } else {
262  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
263  }
264 
265  } // BuildConstantP
266 
267  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
268  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
269  BuildLinearP(RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
270  RCP<realvaluedmultivector_type>& fineCoordinates,
271  RCP<realvaluedmultivector_type>& ghostCoordinates,
272  const int numDimensions, RCP<Matrix>& P) const {
273 
274  // Set debug outputs based on environment variable
275  RCP<Teuchos::FancyOStream> out;
276  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
277  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
278  out->setShowAllFrontMatter(false).setShowProcRank(true);
279  } else {
280  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
281  }
282 
283  *out << "Entering BuildLinearP" << std::endl;
284 
285  // Extract coordinates for interpolation stencil calculations
286  const LO numFineNodes = fineCoordinates->getLocalLength();
287  const LO numGhostNodes = ghostCoordinates->getLocalLength();
288  Array<ArrayRCP<const real_type> > fineCoords(3);
289  Array<ArrayRCP<const real_type> > ghostCoords(3);
290  const real_type realZero = Teuchos::as<real_type>(0.0);
291  ArrayRCP<real_type> fineZero(numFineNodes, realZero);
292  ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
293  for(int dim = 0; dim < 3; ++dim) {
294  if(dim < numDimensions) {
295  fineCoords[dim] = fineCoordinates->getData(dim);
296  ghostCoords[dim] = ghostCoordinates->getData(dim);
297  } else {
298  fineCoords[dim] = fineZero;
299  ghostCoords[dim] = ghostZero;
300  }
301  }
302 
303  *out << "Coordinates extracted from the multivectors!" << std::endl;
304 
305  // Compute 2^numDimensions using bit logic to avoid round-off errors
306  const int numInterpolationPoints = 1 << numDimensions;
307  const int dofsPerNode = A->GetFixedBlockSize();
308 
309  std::vector<size_t> strideInfo(1);
310  strideInfo[0] = dofsPerNode;
311  RCP<const StridedMap> stridedDomainMap =
312  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
313 
314  *out << "The maps of P have been computed" << std::endl;
315 
316  RCP<ParameterList> dummyList = rcp(new ParameterList());
317  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
318  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
319  PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
320 
321  LO interpolationNodeIdx = 0, rowIdx = 0;
322  ArrayView<const LO> colIndices;
323  Array<SC> values;
324  Array<Array<real_type> > coords(numInterpolationPoints + 1);
325  Array<real_type> stencil(numInterpolationPoints);
326  for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
327  if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
328  values.resize(1);
329  values[0] = 1.0;
330  for(LO dof = 0; dof < dofsPerNode; ++dof) {
331  rowIdx = nodeIdx*dofsPerNode + dof;
332  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
333  PCrs->replaceLocalValues(rowIdx, colIndices, values());
334  }
335  } else {
336  // Extract the coordinates associated with the current node
337  // and the neighboring coarse nodes
338  coords[0].resize(3);
339  for(int dim = 0; dim < 3; ++dim) {
340  coords[0][dim] = fineCoords[dim][nodeIdx];
341  }
342  prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
343  for(int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
344  coords[interpolationIdx + 1].resize(3);
345  interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
346  for(int dim = 0; dim < 3; ++dim) {
347  coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
348  }
349  }
350  ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
351  values.resize(numInterpolationPoints);
352  for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
353  values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
354  }
355 
356  // Set values in all the rows corresponding to nodeIdx
357  for(LO dof = 0; dof < dofsPerNode; ++dof) {
358  rowIdx = nodeIdx*dofsPerNode + dof;
359  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
360  PCrs->replaceLocalValues(rowIdx, colIndices, values());
361  }
362  }
363  }
364 
365  *out << "The calculation of the interpolation stencils has completed." << std::endl;
366 
367  PCrs->fillComplete();
368 
369  *out << "All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
370 
371  // set StridingInformation of P
372  if (A->IsView("stridedMaps") == true) {
373  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
374  } else {
375  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
376  }
377 
378  } // BuildLinearP
379 
380 
381  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
382  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
383  ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
384  const Array<Array<real_type> > coord,
385  Array<real_type>& stencil) const {
386 
387  // 7 8 Find xi, eta and zeta such that
388  // x---------x
389  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
390  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
391  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
392  // | | *P | |
393  // | x------|--x We can do this with a Newton solver:
394  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
395  // |/ |/ We compute the Jacobian and iterate until convergence...
396  // z y x---------x
397  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
398  // |/ give us the weights for the interpolation stencil!
399  // o---x
400  //
401 
402  Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
403  Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
404  Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
405  Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
406  Teuchos::SerialDenseSolver<LO,real_type> problem;
407  int iter = 0, max_iter = 5;
408  real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
409  paramCoords.size(numDimensions);
410 
411  while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
412  ++iter;
413  norm2 = 0.0;
414  solutionDirection.size(numDimensions);
415  residual.size(numDimensions);
416  Jacobian = 0.0;
417 
418  // Compute Jacobian and Residual
419  GetInterpolationFunctions(numDimensions, paramCoords, functions);
420  for(LO i = 0; i < numDimensions; ++i) {
421  residual(i) = coord[0][i]; // Add coordinates from point of interest
422  for(LO k = 0; k < numInterpolationPoints; ++k) {
423  residual(i) -= functions[0][k]*coord[k+1][i]; //Remove contribution from all coarse points
424  }
425  if(iter == 1) {
426  norm_ref += residual(i)*residual(i);
427  if(i == numDimensions - 1) {
428  norm_ref = std::sqrt(norm_ref);
429  }
430  }
431 
432  for(LO j = 0; j < numDimensions; ++j) {
433  for(LO k = 0; k < numInterpolationPoints; ++k) {
434  Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
435  }
436  }
437  }
438 
439  // Set Jacobian, Vectors and solve problem
440  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
441  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
442  if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(true);}
443  problem.solve();
444 
445  for(LO i = 0; i < numDimensions; ++i) {
446  paramCoords(i) = paramCoords(i) + solutionDirection(i);
447  }
448 
449  // Recompute Residual norm
450  GetInterpolationFunctions(numDimensions, paramCoords, functions);
451  for(LO i = 0; i < numDimensions; ++i) {
452  real_type tmp = coord[0][i];
453  for(LO k = 0; k < numInterpolationPoints; ++k) {
454  tmp -= functions[0][k]*coord[k+1][i];
455  }
456  norm2 += tmp*tmp;
457  tmp = 0.0;
458  }
459  norm2 = std::sqrt(norm2);
460  }
461 
462  // Load the interpolation values onto the stencil.
463  for(LO i = 0; i < numInterpolationPoints; ++i) {
464  stencil[i] = functions[0][i];
465  }
466 
467  } // End ComputeLinearInterpolationStencil
468 
469  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
470  void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
471  GetInterpolationFunctions(const LO numDimensions,
472  const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
473  real_type functions[4][8]) const {
474  real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
475  if(numDimensions == 1) {
476  xi = parametricCoordinates[0];
477  denominator = 2.0;
478  } else if(numDimensions == 2) {
479  xi = parametricCoordinates[0];
480  eta = parametricCoordinates[1];
481  denominator = 4.0;
482  } else if(numDimensions == 3) {
483  xi = parametricCoordinates[0];
484  eta = parametricCoordinates[1];
485  zeta = parametricCoordinates[2];
486  denominator = 8.0;
487  }
488 
489  functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
490  functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
491  functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
492  functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
493  functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
494  functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
495  functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
496  functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
497 
498  functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
499  functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
500  functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
501  functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
502  functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
503  functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
504  functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
505  functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
506 
507  functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
508  functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
509  functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
510  functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
511  functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
512  functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
513  functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
514  functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
515 
516  functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
517  functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
518  functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
519  functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
520  functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
521  functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
522  functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
523  functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
524 
525  } // End GetInterpolationFunctions
526 
527 } // namespace MueLu
528 
529 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Teuchos::basic_FancyOStream::setShowProcRank
basic_FancyOStream & setShowProcRank(const bool showProcRank)
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::basic_oblackholestream
LO
LocalOrdinal LO
MueLu_GeometricInterpolationPFactory_decl.hpp
MueLu::GeometricInterpolationPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:471
Xpetra_CrsGraph.hpp
MueLu_MasterList.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::GeometricInterpolationPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:115
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::GeometricInterpolationPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:383
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:261
Teuchos::SerialDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
MueLu::GeometricInterpolationPFactory::BuildLinearP
void BuildLinearP(RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, RCP< Matrix > &P) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:269
Teuchos::basic_FancyOStream::setShowAllFrontMatter
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::GeometricInterpolationPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:95
iter
TransListIter iter
Teuchos::SerialDenseVector::size
int size(OrdinalType length_in)
Teuchos::SerialDenseSolver::setVectors
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
MueLu::GeometricInterpolationPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:62
MueLu_Monitor.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::GeometricInterpolationPFactory::BuildConstantP
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:231
MueLu::GeometricInterpolationPFactory::real_type
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
Definition: MueLu_GeometricInterpolationPFactory_decl.hpp:70
Xpetra_CrsMatrixUtils.hpp
MueLu::GeometricInterpolationPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:121
Teuchos::SerialDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)