MueLu  Version of the Day
MueLu_GeometricInterpolationPFactory_kokkos_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_KOKKOS_DEF_HPP
47 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_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_IndexManager_kokkos.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("interp: build coarse coordinates");
67 #undef SET_VALID_ENTRY
68 
69  // general variables needed in GeometricInterpolationPFactory_kokkos
70  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
71  "Generating factory of the matrix A");
72  validParamList->set<RCP<const FactoryBase> >("prolongatorGraph", Teuchos::null,
73  "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
74  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
75  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
76  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
77  "Fine level nullspace used to construct the coarse level nullspace.");
78  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
79  "Number of spacial dimensions in the problem.");
80  validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
81  "Number of nodes per spatial dimension on the coarse grid.");
82  validParamList->set<RCP<const FactoryBase> >("indexManager", Teuchos::null,
83  "The index manager associated with the local mesh.");
84  validParamList->set<RCP<const FactoryBase> >("structuredInterpolationOrder", Teuchos::null,
85  "Interpolation order for constructing the prolongator.");
86 
87  return validParamList;
88  }
89 
90  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
92  DeclareInput(Level& fineLevel, Level& coarseLevel) const {
93  const ParameterList& pL = GetParameterList();
94 
95  Input(fineLevel, "A");
96  Input(fineLevel, "Nullspace");
97  Input(fineLevel, "numDimensions");
98  Input(fineLevel, "prolongatorGraph");
99  Input(fineLevel, "lCoarseNodesPerDim");
100  Input(fineLevel, "structuredInterpolationOrder");
101 
102  if( pL.get<bool>("interp: build coarse coordinates") ||
103  Get<int>(fineLevel, "structuredInterpolationOrder") == 1) {
104  Input(fineLevel, "Coordinates");
105  Input(fineLevel, "indexManager");
106  }
107 
108  }
109 
110  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
112  Build(Level& fineLevel, Level &coarseLevel) const {
113  return BuildP(fineLevel, coarseLevel);
114  }
115 
116  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
118  BuildP(Level &fineLevel, Level &coarseLevel) const {
119  FactoryMonitor m(*this, "BuildP", coarseLevel);
120 
121  // Set debug outputs based on environment variable
122  RCP<Teuchos::FancyOStream> out;
123  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
124  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
125  out->setShowAllFrontMatter(false).setShowProcRank(true);
126  } else {
127  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
128  }
129 
130  *out << "Starting GeometricInterpolationPFactory_kokkos::BuildP." << std::endl;
131 
132  // Get inputs from the parameter list
133  const ParameterList& pL = GetParameterList();
134  const bool buildCoarseCoordinates = pL.get<bool>("interp: build coarse coordinates");
135  const int interpolationOrder = Get<int>(fineLevel, "structuredInterpolationOrder");
136  const int numDimensions = Get<int>(fineLevel, "numDimensions");
137 
138  // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
139  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
140  RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
141  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
142  RCP<Matrix> P;
143 
144  // Check if we need to build coarse coordinates as they are used if we construct
145  // a linear interpolation prolongator
146  if(buildCoarseCoordinates || (interpolationOrder == 1)) {
147  SubFactoryMonitor sfm(*this, "BuildCoordinates", coarseLevel);
148 
149  // Extract data from fine level
150  RCP<IndexManager_kokkos> geoData = Get<RCP<IndexManager_kokkos> >(fineLevel, "indexManager");
151  fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
152 
153  // Build coarse coordinates map/multivector
154  RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
155  Teuchos::OrdinalTraits<GO>::invalid(),
156  geoData->getNumCoarseNodes(),
157  fineCoordinates->getMap()->getIndexBase(),
158  fineCoordinates->getMap()->getComm());
159  coarseCoordinates = Xpetra::MultiVectorFactory<real_type,LO,GO,Node>::
160  Build(coarseCoordsMap, fineCoordinates->getNumVectors());
161 
162  // Construct and launch functor to fill coarse coordinates values
163  // function should take a const view really
164  coarseCoordinatesBuilderFunctor myCoarseCoordinatesBuilder(geoData,
165  fineCoordinates-> getDeviceLocalView(Xpetra::Access::ReadWrite),
166  coarseCoordinates->getDeviceLocalView(Xpetra::Access::OverwriteAll));
167  Kokkos::parallel_for("GeometricInterpolation: build coarse coordinates",
168  Kokkos::RangePolicy<execution_space>(0, geoData->getNumCoarseNodes()),
169  myCoarseCoordinatesBuilder);
170 
171  Set(coarseLevel, "Coordinates", coarseCoordinates);
172  }
173 
174  *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
175 
176  if(interpolationOrder == 0) {
177  SubFactoryMonitor sfm(*this, "BuildConstantP", coarseLevel);
178  // Compute the prolongator using piece-wise constant interpolation
179  BuildConstantP(P, prolongatorGraph, A);
180  } else if(interpolationOrder == 1) {
181  // Compute the prolongator using piece-wise linear interpolation
182  // First get all the required coordinates to compute the local part of P
183  RCP<realvaluedmultivector_type> ghostCoordinates
184  = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(prolongatorGraph->getColMap(),
185  fineCoordinates->getNumVectors());
186  RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
187  prolongatorGraph->getColMap());
188  ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
189 
190  SubFactoryMonitor sfm(*this, "BuildLinearP", coarseLevel);
191  BuildLinearP(A, prolongatorGraph, fineCoordinates, ghostCoordinates, numDimensions, P);
192  }
193 
194  *out << "The prolongator matrix has been built." << std::endl;
195 
196  {
197  SubFactoryMonitor sfm(*this, "BuildNullspace", coarseLevel);
198  // Build the coarse nullspace
199  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
200  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
201  fineNullspace->getNumVectors());
202  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
203  Teuchos::ScalarTraits<SC>::zero());
204  Set(coarseLevel, "Nullspace", coarseNullspace);
205  }
206 
207  *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
208 
209  Array<LO> lNodesPerDir = Get<Array<LO> >(fineLevel, "lCoarseNodesPerDim");
210  Set(coarseLevel, "numDimensions", numDimensions);
211  Set(coarseLevel, "lNodesPerDim", lNodesPerDir);
212  Set(coarseLevel, "P", P);
213 
214  *out << "GeometricInterpolationPFactory_kokkos::BuildP has completed." << std::endl;
215 
216  } // BuildP
217 
218  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
220  BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
221 
222  // Set debug outputs based on environment variable
223  RCP<Teuchos::FancyOStream> out;
224  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
225  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
226  out->setShowAllFrontMatter(false).setShowProcRank(true);
227  } else {
228  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
229  }
230 
231  *out << "BuildConstantP" << std::endl;
232 
233  std::vector<size_t> strideInfo(1);
234  strideInfo[0] = A->GetFixedBlockSize();
235  RCP<const StridedMap> stridedDomainMap =
236  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
237 
238  *out << "Call prolongator constructor" << std::endl;
239 
240  // Create the prolongator matrix and its associated objects
241  RCP<ParameterList> dummyList = rcp(new ParameterList());
242  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
243  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
244  PCrs->setAllToScalar(1.0);
245  PCrs->fillComplete();
246 
247  // set StridingInformation of P
248  if (A->IsView("stridedMaps") == true) {
249  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
250  } else {
251  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
252  }
253 
254  } // BuildConstantP
255 
256  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
258  BuildLinearP(RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
259  RCP<realvaluedmultivector_type>& fineCoordinates,
260  RCP<realvaluedmultivector_type>& ghostCoordinates,
261  const int numDimensions, RCP<Matrix>& P) const {
262 
263  // Set debug outputs based on environment variable
264  RCP<Teuchos::FancyOStream> out;
265  if(const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
266  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
267  out->setShowAllFrontMatter(false).setShowProcRank(true);
268  } else {
269  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
270  }
271 
272  *out << "Entering BuildLinearP" << std::endl;
273 
274  // Extract coordinates for interpolation stencil calculations
275  const LO numFineNodes = fineCoordinates->getLocalLength();
276  const LO numGhostNodes = ghostCoordinates->getLocalLength();
277  Array<ArrayRCP<const real_type> > fineCoords(3);
278  Array<ArrayRCP<const real_type> > ghostCoords(3);
279  const real_type realZero = Teuchos::as<real_type>(0.0);
280  ArrayRCP<real_type> fineZero(numFineNodes, realZero);
281  ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
282  for(int dim = 0; dim < 3; ++dim) {
283  if(dim < numDimensions) {
284  fineCoords[dim] = fineCoordinates->getData(dim);
285  ghostCoords[dim] = ghostCoordinates->getData(dim);
286  } else {
287  fineCoords[dim] = fineZero;
288  ghostCoords[dim] = ghostZero;
289  }
290  }
291 
292  *out << "Coordinates extracted from the multivectors!" << std::endl;
293 
294  // Compute 2^numDimensions using bit logic to avoid round-off errors
295  const int numInterpolationPoints = 1 << numDimensions;
296  const int dofsPerNode = A->GetFixedBlockSize();
297 
298  std::vector<size_t> strideInfo(1);
299  strideInfo[0] = dofsPerNode;
300  RCP<const StridedMap> stridedDomainMap =
301  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
302 
303  *out << "The maps of P have been computed" << std::endl;
304 
305  RCP<ParameterList> dummyList = rcp(new ParameterList());
306  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
307  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
308  PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
309 
310  LO interpolationNodeIdx = 0, rowIdx = 0;
311  ArrayView<const LO> colIndices;
312  Array<SC> values;
313  Array<Array<real_type> > coords(numInterpolationPoints + 1);
314  Array<real_type> stencil(numInterpolationPoints);
315  for(LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
316  if(PCrs->getNumEntriesInLocalRow(nodeIdx*dofsPerNode) == 1) {
317  values.resize(1);
318  values[0] = 1.0;
319  for(LO dof = 0; dof < dofsPerNode; ++dof) {
320  rowIdx = nodeIdx*dofsPerNode + dof;
321  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
322  PCrs->replaceLocalValues(rowIdx, colIndices, values());
323  }
324  } else {
325  // Extract the coordinates associated with the current node
326  // and the neighboring coarse nodes
327  coords[0].resize(3);
328  for(int dim = 0; dim < 3; ++dim) {
329  coords[0][dim] = fineCoords[dim][nodeIdx];
330  }
331  prolongatorGraph->getLocalRowView(nodeIdx*dofsPerNode, colIndices);
332  for(int interpolationIdx=0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
333  coords[interpolationIdx + 1].resize(3);
334  interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
335  for(int dim = 0; dim < 3; ++dim) {
336  coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
337  }
338  }
339  ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
340  values.resize(numInterpolationPoints);
341  for(LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
342  values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
343  }
344 
345  // Set values in all the rows corresponding to nodeIdx
346  for(LO dof = 0; dof < dofsPerNode; ++dof) {
347  rowIdx = nodeIdx*dofsPerNode + dof;
348  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
349  PCrs->replaceLocalValues(rowIdx, colIndices, values());
350  }
351  }
352  }
353 
354  *out << "The calculation of the interpolation stencils has completed." << std::endl;
355 
356  PCrs->fillComplete();
357 
358  *out << "All values in P have been set and expertStaticFillComplete has been performed." << std::endl;
359 
360  // set StridingInformation of P
361  if (A->IsView("stridedMaps") == true) {
362  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
363  } else {
364  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
365  }
366 
367  } // BuildLinearP
368 
369 
370  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
372  ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
373  const Array<Array<real_type> > coord,
374  Array<real_type>& stencil) const {
375 
376  // 7 8 Find xi, eta and zeta such that
377  // x---------x
378  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
379  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
380  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
381  // | | *P | |
382  // | x------|--x We can do this with a Newton solver:
383  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
384  // |/ |/ We compute the Jacobian and iterate until convergence...
385  // z y x---------x
386  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
387  // |/ give us the weights for the interpolation stencil!
388  // o---x
389  //
390 
391  Teuchos::SerialDenseMatrix<LO,real_type> Jacobian(numDimensions, numDimensions);
392  Teuchos::SerialDenseVector<LO,real_type> residual(numDimensions);
393  Teuchos::SerialDenseVector<LO,real_type> solutionDirection(numDimensions);
394  Teuchos::SerialDenseVector<LO,real_type> paramCoords(numDimensions);
395  Teuchos::SerialDenseSolver<LO,real_type> problem;
396  int iter = 0, max_iter = 5;
397  real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
398  paramCoords.size(numDimensions);
399 
400  while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
401  ++iter;
402  norm2 = 0.0;
403  solutionDirection.size(numDimensions);
404  residual.size(numDimensions);
405  Jacobian = 0.0;
406 
407  // Compute Jacobian and Residual
408  GetInterpolationFunctions(numDimensions, paramCoords, functions);
409  for(LO i = 0; i < numDimensions; ++i) {
410  residual(i) = coord[0][i]; // Add coordinates from point of interest
411  for(LO k = 0; k < numInterpolationPoints; ++k) {
412  residual(i) -= functions[0][k]*coord[k+1][i]; //Remove contribution from all coarse points
413  }
414  if(iter == 1) {
415  norm_ref += residual(i)*residual(i);
416  if(i == numDimensions - 1) {
417  norm_ref = std::sqrt(norm_ref);
418  }
419  }
420 
421  for(LO j = 0; j < numDimensions; ++j) {
422  for(LO k = 0; k < numInterpolationPoints; ++k) {
423  Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
424  }
425  }
426  }
427 
428  // Set Jacobian, Vectors and solve problem
429  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
430  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
431  if(problem.shouldEquilibrate()) {problem.factorWithEquilibration(true);}
432  problem.solve();
433 
434  for(LO i = 0; i < numDimensions; ++i) {
435  paramCoords(i) = paramCoords(i) + solutionDirection(i);
436  }
437 
438  // Recompute Residual norm
439  GetInterpolationFunctions(numDimensions, paramCoords, functions);
440  for(LO i = 0; i < numDimensions; ++i) {
441  real_type tmp = coord[0][i];
442  for(LO k = 0; k < numInterpolationPoints; ++k) {
443  tmp -= functions[0][k]*coord[k+1][i];
444  }
445  norm2 += tmp*tmp;
446  tmp = 0.0;
447  }
448  norm2 = std::sqrt(norm2);
449  }
450 
451  // Load the interpolation values onto the stencil.
452  for(LO i = 0; i < numInterpolationPoints; ++i) {
453  stencil[i] = functions[0][i];
454  }
455 
456  } // End ComputeLinearInterpolationStencil
457 
458  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
460  GetInterpolationFunctions(const LO numDimensions,
461  const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
462  real_type functions[4][8]) const {
463  real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
464  if(numDimensions == 1) {
465  xi = parametricCoordinates[0];
466  denominator = 2.0;
467  } else if(numDimensions == 2) {
468  xi = parametricCoordinates[0];
469  eta = parametricCoordinates[1];
470  denominator = 4.0;
471  } else if(numDimensions == 3) {
472  xi = parametricCoordinates[0];
473  eta = parametricCoordinates[1];
474  zeta = parametricCoordinates[2];
475  denominator = 8.0;
476  }
477 
478  functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
479  functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
480  functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
481  functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
482  functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
483  functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
484  functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
485  functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
486 
487  functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
488  functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
489  functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
490  functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
491  functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
492  functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
493  functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
494  functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
495 
496  functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
497  functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
498  functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
499  functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
500  functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
501  functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
502  functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
503  functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
504 
505  functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
506  functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
507  functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
508  functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
509  functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
510  functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
511  functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
512  functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
513 
514  } // End GetInterpolationFunctions
515 
516  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
519  coord_view_type fineCoordView,
520  coord_view_type coarseCoordView)
521  : geoData_(*geoData), fineCoordView_(fineCoordView), coarseCoordView_(coarseCoordView) {
522 
523  } // coarseCoordinatesBuilderFunctor()
524 
525  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
526  KOKKOS_INLINE_FUNCTION
528  coarseCoordinatesBuilderFunctor::operator() (const LO coarseNodeIdx) const {
529 
530  LO fineNodeIdx;
531  LO nodeCoarseTuple[3] = {0, 0, 0};
532  LO nodeFineTuple[3] = {0, 0, 0};
533  auto coarseningRate = geoData_.getCoarseningRates();
534  auto fineNodesPerDir = geoData_.getLocalFineNodesPerDir();
535  auto coarseNodesPerDir = geoData_.getCoarseNodesPerDir();
536  geoData_.getCoarseLID2CoarseTuple(coarseNodeIdx, nodeCoarseTuple);
537  for(int dim = 0; dim < 3; ++dim) {
538  if(nodeCoarseTuple[dim] == coarseNodesPerDir(dim) - 1) {
539  nodeFineTuple[dim] = fineNodesPerDir(dim) - 1;
540  } else {
541  nodeFineTuple[dim] = nodeCoarseTuple[dim]*coarseningRate(dim);
542  }
543  }
544 
545  fineNodeIdx = nodeFineTuple[2]*fineNodesPerDir(1)*fineNodesPerDir(0)
546  + nodeFineTuple[1]*fineNodesPerDir(0) + nodeFineTuple[0];
547 
548  for(int dim = 0; dim < fineCoordView_.extent_int(1); ++dim) {
549  coarseCoordView_(coarseNodeIdx, dim) = fineCoordView_(fineNodeIdx, dim);
550  }
551  }
552 
553 } // namespace MueLu
554 
555 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_KOKKOS_DEF_HPP
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
coarseCoordinatesBuilderFunctor(RCP< IndexManager_kokkos > geoData, coord_view_type fineCoordView, coord_view_type coarseCoordView)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Namespace for MueLu classes and methods.
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
void BuildLinearP(RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, RCP< Matrix > &P) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
typename Kokkos::View< impl_scalar_type **, Kokkos::LayoutLeft, device_type > coord_view_type
#define SET_VALID_ENTRY(name)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const