MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_GeneralGeometricPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
11 #define MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
12 
13 #include <stdlib.h>
14 #include <iomanip>
15 
16 // #include <Teuchos_LAPACK.hpp>
17 #include <Teuchos_SerialDenseMatrix.hpp>
18 #include <Teuchos_SerialDenseVector.hpp>
19 #include <Teuchos_SerialDenseSolver.hpp>
20 
21 #include <Xpetra_CrsMatrixWrap.hpp>
22 #include <Xpetra_ImportFactory.hpp>
23 #include <Xpetra_Matrix.hpp>
24 #include <Xpetra_MapFactory.hpp>
25 #include <Xpetra_MultiVectorFactory.hpp>
26 #include <Xpetra_VectorFactory.hpp>
27 
28 #include <Xpetra_IO.hpp>
29 
30 #include "MueLu_GeneralGeometricPFactory_decl.hpp"
31 
32 #include "MueLu_Monitor.hpp"
33 
34 namespace MueLu {
35 
36 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
37 RCP<const ParameterList> GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
38  RCP<ParameterList> validParamList = rcp(new ParameterList());
39 
40  // Coarsen can come in two forms, either a single char that will be interpreted as an integer
41  // which is used as the coarsening rate in every spatial dimentions, or it can be a longer
42  // string that will then be interpreted as an array of integers.
43  // By default coarsen is set as "{2}", hence a coarsening rate of 2 in every spatial dimension
44  // is the default setting!
45  validParamList->set<std::string>("Coarsen", "{2}",
46  "Coarsening rate in all spatial dimensions");
47  validParamList->set<int>("order", 1,
48  "Order of the interpolation scheme used");
49  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
50  "Generating factory of the matrix A");
51  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
52  "Generating factory of the nullspace");
53  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
54  "Generating factory for coorindates");
55  validParamList->set<RCP<const FactoryBase> >("gNodesPerDim", Teuchos::null,
56  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
57  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
58  "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
59  validParamList->set<std::string>("meshLayout", "Global Lexicographic",
60  "Type of mesh ordering");
61  validParamList->set<RCP<const FactoryBase> >("meshData", Teuchos::null,
62  "Mesh ordering associated data");
63 
64  return validParamList;
65 }
66 
67 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
68 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
69  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
70  Input(fineLevel, "A");
71  Input(fineLevel, "Nullspace");
72  Input(fineLevel, "Coordinates");
73  // Request the global number of nodes per dimensions
74  if (fineLevel.GetLevelID() == 0) {
75  if (fineLevel.IsAvailable("gNodesPerDim", NoFactory::get())) {
76  fineLevel.DeclareInput("gNodesPerDim", NoFactory::get(), this);
77  } else {
78  TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("gNodesPerDim", NoFactory::get()),
79  Exceptions::RuntimeError,
80  "gNodesPerDim was not provided by the user on level0!");
81  }
82  } else {
83  Input(fineLevel, "gNodesPerDim");
84  }
85 
86  // Request the local number of nodes per dimensions
87  if (fineLevel.GetLevelID() == 0) {
88  if (fineLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
89  fineLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
90  } else {
91  TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
92  Exceptions::RuntimeError,
93  "lNodesPerDim was not provided by the user on level0!");
94  }
95  } else {
96  Input(fineLevel, "lNodesPerDim");
97  }
98 }
99 
100 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
101 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel,
102  Level& coarseLevel) const {
103  return BuildP(fineLevel, coarseLevel);
104 }
105 
106 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
107 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& fineLevel,
108  Level& coarseLevel) const {
109  FactoryMonitor m(*this, "Build", coarseLevel);
110 
111  // Obtain general variables
112  using xdMV = Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>;
113  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
114  RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel, "Nullspace");
115  RCP<xdMV> fineCoords = Get<RCP<xdMV> >(fineLevel, "Coordinates");
116  RCP<xdMV> coarseCoords;
117 
118  // Get user-provided coarsening rate parameter (constant over all levels)
119  const ParameterList& pL = GetParameterList();
120 
121  // collect general input data
122  const LO blkSize = A->GetFixedBlockSize();
123  RCP<const Map> rowMap = A->getRowMap();
124  RCP<GeometricData> myGeometry = rcp(new GeometricData{});
125 
126  // Load the mesh layout type and the associated mesh data
127  myGeometry->meshLayout = pL.get<std::string>("meshLayout");
128  if (fineLevel.GetLevelID() == 0) {
129  if (myGeometry->meshLayout == "Local Lexicographic") {
130  Array<GO> meshData;
131  meshData = fineLevel.Get<Array<GO> >("meshData", NoFactory::get());
132  TEUCHOS_TEST_FOR_EXCEPTION(meshData.empty() == true, Exceptions::RuntimeError,
133  "The meshData array is empty, somehow the input for geometric"
134  " multigrid are not captured correctly.");
135  myGeometry->meshData.resize(rowMap->getComm()->getSize());
136  for (int i = 0; i < rowMap->getComm()->getSize(); ++i) {
137  myGeometry->meshData[i].resize(10);
138  for (int j = 0; j < 10; ++j) {
139  myGeometry->meshData[i][j] = meshData[10 * i + j];
140  }
141  }
142  }
143  }
144 
145  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords == Teuchos::null, Exceptions::RuntimeError,
146  "Coordinates cannot be accessed from fine level!");
147  myGeometry->numDimensions = fineCoords->getNumVectors();
148 
149  // Get the number of points in each direction
150  if (fineLevel.GetLevelID() == 0) {
151  myGeometry->gFineNodesPerDir = fineLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
152  myGeometry->lFineNodesPerDir = fineLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
153  } else {
154  // Loading global number of nodes per diretions
155  myGeometry->gFineNodesPerDir = Get<Array<GO> >(fineLevel, "gNodesPerDim");
156 
157  // Loading local number of nodes per diretions
158  myGeometry->lFineNodesPerDir = Get<Array<LO> >(fineLevel, "lNodesPerDim");
159  }
160  myGeometry->gNumFineNodes10 = myGeometry->gFineNodesPerDir[1] * myGeometry->gFineNodesPerDir[0];
161  myGeometry->gNumFineNodes = myGeometry->gFineNodesPerDir[2] * myGeometry->gNumFineNodes10;
162  myGeometry->lNumFineNodes10 = myGeometry->lFineNodesPerDir[1] * myGeometry->lFineNodesPerDir[0];
163  myGeometry->lNumFineNodes = myGeometry->lFineNodesPerDir[2] * myGeometry->lNumFineNodes10;
164 
165  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getLocalLength() != static_cast<size_t>(myGeometry->lNumFineNodes),
166  Exceptions::RuntimeError,
167  "The local number of elements in Coordinates is not equal to the"
168  " number of nodes given by: lNodesPerDim!");
169  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords->getGlobalLength() != static_cast<size_t>(myGeometry->gNumFineNodes),
170  Exceptions::RuntimeError,
171  "The global number of elements in Coordinates is not equal to the"
172  " number of nodes given by: gNodesPerDim!");
173 
174  // Get the coarsening rate
175  std::string coarsenRate = pL.get<std::string>("Coarsen");
176  Teuchos::Array<LO> crates;
177  try {
178  crates = Teuchos::fromStringToArray<LO>(coarsenRate);
179  } catch (const Teuchos::InvalidArrayStringRepresentation& e) {
180  GetOStream(Errors, -1) << " *** Coarsen must be a string convertible into an array! *** "
181  << std::endl;
182  throw e;
183  }
184  TEUCHOS_TEST_FOR_EXCEPTION((crates.size() > 1) && (crates.size() < myGeometry->numDimensions),
185  Exceptions::RuntimeError,
186  "Coarsen must have at least as many components as the number of"
187  " spatial dimensions in the problem.");
188 
189  for (LO i = 0; i < 3; ++i) {
190  if (i < myGeometry->numDimensions) {
191  if (crates.size() == 1) {
192  myGeometry->coarseRate[i] = crates[0];
193  } else if (crates.size() == myGeometry->numDimensions) {
194  myGeometry->coarseRate[i] = crates[i];
195  }
196  } else {
197  myGeometry->coarseRate[i] = 1;
198  }
199  }
200 
201  int interpolationOrder = pL.get<int>("order");
202  TEUCHOS_TEST_FOR_EXCEPTION((interpolationOrder < 0) || (interpolationOrder > 1),
203  Exceptions::RuntimeError,
204  "The interpolation order can only be set to 0 or 1.");
205 
206  // Get the axis permutation from Global axis to Local axis
207  Array<LO> mapDirG2L(3), mapDirL2G(3);
208  for (LO i = 0; i < myGeometry->numDimensions; ++i) {
209  mapDirG2L[i] = i;
210  }
211  for (LO i = 0; i < myGeometry->numDimensions; ++i) {
212  TEUCHOS_TEST_FOR_EXCEPTION(mapDirG2L[i] > myGeometry->numDimensions,
213  Exceptions::RuntimeError,
214  "axis permutation values must all be less than"
215  " the number of spatial dimensions.");
216  mapDirL2G[mapDirG2L[i]] = i;
217  }
218  RCP<const Map> fineCoordsMap = fineCoords->getMap();
219 
220  // This struct stores PIDs, LIDs and GIDs on the fine mesh and GIDs on the coarse mesh.
221  RCP<NodesIDs> ghostedCoarseNodes = rcp(new NodesIDs{});
222  Array<Array<GO> > lCoarseNodesGIDs(2);
223  if ((fineLevel.GetLevelID() == 0) && (myGeometry->meshLayout != "Global Lexicographic")) {
224  MeshLayoutInterface(interpolationOrder, blkSize, fineCoordsMap, myGeometry,
225  ghostedCoarseNodes, lCoarseNodesGIDs);
226  } else {
227  // This function expects perfect global lexicographic ordering of nodes and will not process
228  // data correctly otherwise. These restrictions allow for the simplest and most efficient
229  // processing of the levels (hopefully at least).
230  GetCoarsePoints(interpolationOrder, blkSize, fineCoordsMap, myGeometry, ghostedCoarseNodes,
231  lCoarseNodesGIDs);
232  }
233 
234  // All that is left to do is loop over NCpts and:
235  // - extract coarse points coordiate for coarseCoords
236  // - get coordinates for current stencil computation
237  // - compute current stencil
238  // - compute row and column indices for stencil entries
239  RCP<const Map> stridedDomainMapP;
240  RCP<Matrix> P;
241  // Fancy formula for the number of non-zero terms
242  // All coarse points are injected, other points are using polynomial interpolation
243  // and have contribution from (interpolationOrder + 1)^numDimensions
244  // Noticebly this leads to 1 when the order is zero, hence fancy MatMatMatMult can be used.
245  GO lnnzP = Teuchos::as<LO>(std::pow(interpolationOrder + 1, myGeometry->numDimensions)) * (myGeometry->lNumFineNodes - myGeometry->lNumCoarseNodes) + myGeometry->lNumCoarseNodes;
246  lnnzP = lnnzP * blkSize;
247  GO gnnzP = Teuchos::as<LO>(std::pow(interpolationOrder + 1, myGeometry->numDimensions)) * (myGeometry->gNumFineNodes - myGeometry->gNumCoarseNodes) + myGeometry->gNumCoarseNodes;
248  gnnzP = gnnzP * blkSize;
249 
250  GetOStream(Runtime1) << "P size = " << blkSize * myGeometry->gNumFineNodes
251  << " x " << blkSize * myGeometry->gNumCoarseNodes << std::endl;
252  GetOStream(Runtime1) << "P Fine grid : " << myGeometry->gFineNodesPerDir[0] << " -- "
253  << myGeometry->gFineNodesPerDir[1] << " -- "
254  << myGeometry->gFineNodesPerDir[2] << std::endl;
255  GetOStream(Runtime1) << "P Coarse grid : " << myGeometry->gCoarseNodesPerDir[0] << " -- "
256  << myGeometry->gCoarseNodesPerDir[1] << " -- "
257  << myGeometry->gCoarseNodesPerDir[2] << std::endl;
258  GetOStream(Runtime1) << "P nnz estimate: " << gnnzP << std::endl;
259 
260  MakeGeneralGeometricP(myGeometry, fineCoords, lnnzP, blkSize, stridedDomainMapP,
261  A, P, coarseCoords, ghostedCoarseNodes, lCoarseNodesGIDs,
262  interpolationOrder);
263 
264  // set StridingInformation of P
265  if (A->IsView("stridedMaps") == true) {
266  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMapP);
267  } else {
268  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMapP);
269  }
270 
271  // store the transfer operator and node coordinates on coarse level
272  Set(coarseLevel, "P", P);
273  Set(coarseLevel, "coarseCoordinates", coarseCoords);
274  Set<Array<GO> >(coarseLevel, "gCoarseNodesPerDim", myGeometry->gCoarseNodesPerDir);
275  Set<Array<LO> >(coarseLevel, "lCoarseNodesPerDim", myGeometry->lCoarseNodesPerDir);
276 
277  // rst: null space might get scaled here ... do we care. We could just inject at the cpoints,
278  // but I don't feel that this is needed.
279  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
280  fineNullspace->getNumVectors());
281  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
282  Teuchos::ScalarTraits<SC>::zero());
283  Set(coarseLevel, "Nullspace", coarseNullspace);
284 }
285 
286 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
287 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
288  MeshLayoutInterface(const int /* interpolationOrder */, const LO /* blkSize */, RCP<const Map> fineCoordsMap,
289  RCP<GeometricData> myGeo, RCP<NodesIDs> ghostedCoarseNodes,
290  Array<Array<GO> >& lCoarseNodesGIDs) const {
291  // The goal here is to produce maps that globally labels the mesh lexicographically.
292  // These maps will replace the current maps of A, the coordinate vector and the nullspace.
293  // Ideally if the user provides the necessary maps then nothing needs to be done, otherwise
294  // it could be advantageous to allow the user to register a re-labeling function. Ultimately
295  // for common ordering schemes, some re-labeling can be directly implemented here.
296 
297  int numRanks = fineCoordsMap->getComm()->getSize();
298  int myRank = fineCoordsMap->getComm()->getRank();
299 
300  myGeo->myBlock = myGeo->meshData[myRank][2];
301  myGeo->startIndices[0] = myGeo->meshData[myRank][3];
302  myGeo->startIndices[1] = myGeo->meshData[myRank][5];
303  myGeo->startIndices[2] = myGeo->meshData[myRank][7];
304  myGeo->startIndices[3] = myGeo->meshData[myRank][4];
305  myGeo->startIndices[4] = myGeo->meshData[myRank][6];
306  myGeo->startIndices[5] = myGeo->meshData[myRank][8];
307  std::sort(myGeo->meshData.begin(), myGeo->meshData.end(),
308  [](const std::vector<GO>& a, const std::vector<GO>& b) -> bool {
309  // The below function sorts ranks by blockID, kmin, jmin and imin
310  if (a[2] < b[2]) {
311  return true;
312  } else if (a[2] == b[2]) {
313  if (a[7] < b[7]) {
314  return true;
315  } else if (a[7] == b[7]) {
316  if (a[5] < b[5]) {
317  return true;
318  } else if (a[5] == b[5]) {
319  if (a[3] < b[3]) {
320  return true;
321  }
322  }
323  }
324  }
325  return false;
326  });
327 
328  myGeo->numBlocks = myGeo->meshData[numRanks - 1][2] + 1;
329  // Find the range of the current block
330  auto myBlockStart = std::lower_bound(myGeo->meshData.begin(), myGeo->meshData.end(),
331  myGeo->myBlock - 1,
332  [](const std::vector<GO>& vec, const GO val) -> bool {
333  return (vec[2] < val) ? true : false;
334  });
335  auto myBlockEnd = std::upper_bound(myGeo->meshData.begin(), myGeo->meshData.end(),
336  myGeo->myBlock,
337  [](const GO val, const std::vector<GO>& vec) -> bool {
338  return (val < vec[2]) ? true : false;
339  });
340  // Assuming that i,j,k and ranges are split in pi, pj and pk processors
341  // we search for these numbers as they will allow us to find quickly the PID of processors
342  // owning ghost nodes.
343  auto myKEnd = std::upper_bound(myBlockStart, myBlockEnd, (*myBlockStart)[3],
344  [](const GO val, const std::vector<GO>& vec) -> bool {
345  return (val < vec[7]) ? true : false;
346  });
347  auto myJEnd = std::upper_bound(myBlockStart, myKEnd, (*myBlockStart)[3],
348  [](const GO val, const std::vector<GO>& vec) -> bool {
349  return (val < vec[5]) ? true : false;
350  });
351  LO pi = std::distance(myBlockStart, myJEnd);
352  LO pj = std::distance(myBlockStart, myKEnd) / pi;
353  LO pk = std::distance(myBlockStart, myBlockEnd) / (pj * pi);
354 
355  // We also look for the index of the local rank in the current block.
356  LO myRankIndex = std::distance(myGeo->meshData.begin(),
357  std::find_if(myBlockStart, myBlockEnd,
358  [myRank](const std::vector<GO>& vec) -> bool {
359  return (vec[0] == myRank) ? true : false;
360  }));
361 
362  for (int dim = 0; dim < 3; ++dim) {
363  if (dim < myGeo->numDimensions) {
364  myGeo->offsets[dim] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
365  myGeo->offsets[dim + 3] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
366  }
367  }
368 
369  // Check if the partition contains nodes on a boundary, if so that boundary (face, line or
370  // point) will not require ghost nodes.
371  for (int dim = 0; dim < 3; ++dim) {
372  if (dim < myGeo->numDimensions && (myGeo->startIndices[dim] % myGeo->coarseRate[dim] != 0 ||
373  myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1)) {
374  myGeo->ghostInterface[2 * dim] = true;
375  }
376  if (dim < myGeo->numDimensions && myGeo->startIndices[dim + 3] != myGeo->gFineNodesPerDir[dim] - 1 && (myGeo->lFineNodesPerDir[dim] == 1 || myGeo->startIndices[dim + 3] % myGeo->coarseRate[dim] != 0)) {
377  myGeo->ghostInterface[2 * dim + 1] = true;
378  }
379  }
380 
381  // Here one element can represent either the degenerate case of one node or the more general
382  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with one
383  // node. This helps generating a 3D space from tensorial products...
384  // A good way to handle this would be to generalize the algorithm to take into account the
385  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
386  // discretization will have a unique node per element. This way 1D discretization can be viewed
387  // as a 3D problem with one 0 degree element in the y direction and one 0 degre element in the z
388  // direction.
389  // !!! Operations below are aftecting both local and global values that have two different !!!
390  // orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G. coarseRate,
391  // endRate and offsets are in the global basis, as well as all the variables starting with a g,
392  // !!! while the variables starting with an l are in the local basis. !!!
393  for (int i = 0; i < 3; ++i) {
394  if (i < myGeo->numDimensions) {
395  // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next
396  // level.
397  myGeo->gCoarseNodesPerDir[i] = (myGeo->gFineNodesPerDir[i] - 1) / myGeo->coarseRate[i];
398  myGeo->endRate[i] = Teuchos::as<LO>((myGeo->gFineNodesPerDir[i] - 1) % myGeo->coarseRate[i]);
399  if (myGeo->endRate[i] == 0) {
400  myGeo->endRate[i] = myGeo->coarseRate[i];
401  ++myGeo->gCoarseNodesPerDir[i];
402  } else {
403  myGeo->gCoarseNodesPerDir[i] += 2;
404  }
405  } else {
406  myGeo->endRate[i] = 1;
407  myGeo->gCoarseNodesPerDir[i] = 1;
408  }
409  }
410 
411  myGeo->gNumCoarseNodes = myGeo->gCoarseNodesPerDir[0] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[2];
412 
413  for (LO i = 0; i < 3; ++i) {
414  if (i < myGeo->numDimensions) {
415  // Check whether the partition includes the "end" of the mesh which means that endRate will
416  // apply. Also make sure that endRate is not 0 which means that the mesh does not require a
417  // particular treatment at the boundaries.
418  if ((myGeo->startIndices[i] + myGeo->lFineNodesPerDir[i]) == myGeo->gFineNodesPerDir[i]) {
419  myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] - myGeo->endRate[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i] + 1;
420  if (myGeo->offsets[i] == 0) {
421  ++myGeo->lCoarseNodesPerDir[i];
422  }
423  } else {
424  myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i];
425  if (myGeo->offsets[i] == 0) {
426  ++myGeo->lCoarseNodesPerDir[i];
427  }
428  }
429  } else {
430  myGeo->lCoarseNodesPerDir[i] = 1;
431  }
432  // This would happen if the rank does not own any nodes but in that case a subcommunicator
433  // should be used so this should really not be a concern.
434  if (myGeo->lFineNodesPerDir[i] < 1) {
435  myGeo->lCoarseNodesPerDir[i] = 0;
436  }
437  }
438 
439  // Assuming linear interpolation, each fine point has contribution from 8 coarse points
440  // and each coarse point value gets injected.
441  // For systems of PDEs we assume that all dofs have the same P operator.
442  myGeo->lNumCoarseNodes = myGeo->lCoarseNodesPerDir[0] * myGeo->lCoarseNodesPerDir[1] * myGeo->lCoarseNodesPerDir[2];
443 
444  // For each direction, determine how many points (including ghosts) are required.
445  for (int dim = 0; dim < 3; ++dim) {
446  // The first branch of this if-statement will be used if the rank contains only one layer
447  // of nodes in direction i, that layer must also coincide with the boundary of the mesh
448  // and coarseRate[i] == endRate[i]...
449  if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1 &&
450  myGeo->startIndices[dim] % myGeo->coarseRate[dim] == 0) {
451  myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim] - 1;
452  } else {
453  myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim];
454  }
455  myGeo->ghostedCoarseNodesPerDir[dim] = myGeo->lCoarseNodesPerDir[dim];
456  // Check whether face *low needs ghost nodes
457  if (myGeo->ghostInterface[2 * dim]) {
458  myGeo->ghostedCoarseNodesPerDir[dim] += 1;
459  }
460  // Check whether face *hi needs ghost nodes
461  if (myGeo->ghostInterface[2 * dim + 1]) {
462  myGeo->ghostedCoarseNodesPerDir[dim] += 1;
463  }
464  }
465  myGeo->lNumGhostedNodes = myGeo->ghostedCoarseNodesPerDir[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0];
466  myGeo->lNumGhostNodes = myGeo->lNumGhostedNodes - myGeo->lNumCoarseNodes;
467  ghostedCoarseNodes->PIDs.resize(myGeo->lNumGhostedNodes);
468  ghostedCoarseNodes->LIDs.resize(myGeo->lNumGhostedNodes);
469  ghostedCoarseNodes->GIDs.resize(myGeo->lNumGhostedNodes);
470  ghostedCoarseNodes->coarseGIDs.resize(myGeo->lNumGhostedNodes);
471  ghostedCoarseNodes->colInds.resize(myGeo->lNumGhostedNodes);
472  lCoarseNodesGIDs[0].resize(myGeo->lNumCoarseNodes);
473  lCoarseNodesGIDs[1].resize(myGeo->lNumCoarseNodes);
474 
475  // Now the tricky part starts, the coarse nodes / ghosted coarse nodes need to be imported.
476  // This requires finding what their GID on the fine mesh is. They need to be ordered
477  // lexicographically to allow for fast sweeps through the mesh.
478 
479  // We loop over all ghosted coarse nodes by increasing global lexicographic order
480  Array<LO> coarseNodeCoarseIndices(3), coarseNodeFineIndices(3);
481  LO currentIndex = -1, countCoarseNodes = 0;
482  for (int k = 0; k < myGeo->ghostedCoarseNodesPerDir[2]; ++k) {
483  for (int j = 0; j < myGeo->ghostedCoarseNodesPerDir[1]; ++j) {
484  for (int i = 0; i < myGeo->ghostedCoarseNodesPerDir[0]; ++i) {
485  currentIndex = k * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + j * myGeo->ghostedCoarseNodesPerDir[0] + i;
486  coarseNodeCoarseIndices[0] = myGeo->startGhostedCoarseNode[0] + i;
487  coarseNodeFineIndices[0] = coarseNodeCoarseIndices[0] * myGeo->coarseRate[0];
488  if (coarseNodeFineIndices[0] > myGeo->gFineNodesPerDir[0] - 1) {
489  coarseNodeFineIndices[0] = myGeo->gFineNodesPerDir[0] - 1;
490  }
491  coarseNodeCoarseIndices[1] = myGeo->startGhostedCoarseNode[1] + j;
492  coarseNodeFineIndices[1] = coarseNodeCoarseIndices[1] * myGeo->coarseRate[1];
493  if (coarseNodeFineIndices[1] > myGeo->gFineNodesPerDir[1] - 1) {
494  coarseNodeFineIndices[1] = myGeo->gFineNodesPerDir[1] - 1;
495  }
496  coarseNodeCoarseIndices[2] = myGeo->startGhostedCoarseNode[2] + k;
497  coarseNodeFineIndices[2] = coarseNodeCoarseIndices[2] * myGeo->coarseRate[2];
498  if (coarseNodeFineIndices[2] > myGeo->gFineNodesPerDir[2] - 1) {
499  coarseNodeFineIndices[2] = myGeo->gFineNodesPerDir[2] - 1;
500  }
501  GO myGID = -1, myCoarseGID = -1;
502  LO myLID = -1, myPID = -1;
503  GetGIDLocalLexicographic(i, j, k, coarseNodeFineIndices, myGeo, myRankIndex, pi, pj, pk,
504  myBlockStart, myBlockEnd, myGID, myPID, myLID);
505  myCoarseGID = coarseNodeCoarseIndices[0] + coarseNodeCoarseIndices[1] * myGeo->gCoarseNodesPerDir[0] + coarseNodeCoarseIndices[2] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[0];
506  ghostedCoarseNodes->PIDs[currentIndex] = myPID;
507  ghostedCoarseNodes->LIDs[currentIndex] = myLID;
508  ghostedCoarseNodes->GIDs[currentIndex] = myGID;
509  ghostedCoarseNodes->coarseGIDs[currentIndex] = myCoarseGID;
510  if (myPID == myRank) {
511  lCoarseNodesGIDs[0][countCoarseNodes] = myCoarseGID;
512  lCoarseNodesGIDs[1][countCoarseNodes] = myGID;
513  ++countCoarseNodes;
514  }
515  }
516  }
517  }
518 
519 } // End MeshLayoutInterface
520 
521 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
522 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
523  GetCoarsePoints(const int /* interpolationOrder */, const LO /* blkSize */, RCP<const Map> fineCoordsMap,
524  RCP<GeometricData> myGeo, RCP<NodesIDs> ghostedCoarseNodes,
525  Array<Array<GO> >& lCoarseNodesGIDs) const {
526  // Assuming perfect global lexicographic ordering of the mesh, produce two arrays:
527  // 1) lGhostNodesIDs that stores PID, LID, GID and coarseGID associated with the coarse nodes
528  // need to compute the local part of the prolongator.
529  // 2) lCoarseNodesGIDs that stores the GIDs associated with the local nodes needed to create
530  // the map of the MultiVector of coarse node coordinates.
531 
532  {
533  GO tmp = 0;
534  myGeo->startIndices[2] = fineCoordsMap->getMinGlobalIndex() / (myGeo->gFineNodesPerDir[1] * myGeo->gFineNodesPerDir[0]);
535  tmp = fineCoordsMap->getMinGlobalIndex() % (myGeo->gFineNodesPerDir[1] * myGeo->gFineNodesPerDir[0]);
536  myGeo->startIndices[1] = tmp / myGeo->gFineNodesPerDir[0];
537  myGeo->startIndices[0] = tmp % myGeo->gFineNodesPerDir[0];
538  } // End of scope for tmp
539  for (int dim = 0; dim < 3; ++dim) {
540  myGeo->startIndices[dim + 3] = myGeo->startIndices[dim] + myGeo->lFineNodesPerDir[dim] - 1;
541  }
542 
543  for (int dim = 0; dim < 3; ++dim) {
544  if (dim < myGeo->numDimensions) {
545  myGeo->offsets[dim] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
546  myGeo->offsets[dim + 3] = Teuchos::as<LO>(myGeo->startIndices[dim]) % myGeo->coarseRate[dim];
547  }
548  }
549 
550  // Check if the partition contains nodes on a boundary, if so that boundary (face, line or
551  // point) will not require ghost nodes, unless there is only one node in that direction.
552  for (int dim = 0; dim < 3; ++dim) {
553  if (dim < myGeo->numDimensions && (myGeo->startIndices[dim] % myGeo->coarseRate[dim] != 0 ||
554  myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1)) {
555  myGeo->ghostInterface[2 * dim] = true;
556  }
557  if (dim < myGeo->numDimensions && myGeo->startIndices[dim + 3] != myGeo->gFineNodesPerDir[dim] - 1 && (myGeo->lFineNodesPerDir[dim] == 1 || myGeo->startIndices[dim + 3] % myGeo->coarseRate[dim] != 0)) {
558  myGeo->ghostInterface[2 * dim + 1] = true;
559  }
560  }
561 
562  // Here one element can represent either the degenerate case of one node or the more general
563  // case of two nodes, i.e. x---x is a 1D element with two nodes and x is a 1D element with one
564  // node. This helps generating a 3D space from tensorial products...
565  // A good way to handle this would be to generalize the algorithm to take into account the
566  // discretization order used in each direction, at least in the FEM sense, since a 0 degree
567  // discretization will have a unique node per element. This way 1D discretization can be viewed
568  // as a 3D problem with one 0 degree element in the y direction and one 0 degre element in the z
569  // direction.
570  // !!! Operations below are aftecting both local and global values that have two different !!!
571  // orientations. Orientations can be interchanged using mapDirG2L and mapDirL2G. coarseRate,
572  // endRate and offsets are in the global basis, as well as all the variables starting with a g,
573  // !!! while the variables starting with an l are in the local basis. !!!
574  for (int i = 0; i < 3; ++i) {
575  if (i < myGeo->numDimensions) {
576  // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next
577  // level.
578  myGeo->gCoarseNodesPerDir[i] = (myGeo->gFineNodesPerDir[i] - 1) / myGeo->coarseRate[i];
579  myGeo->endRate[i] = Teuchos::as<LO>((myGeo->gFineNodesPerDir[i] - 1) % myGeo->coarseRate[i]);
580  if (myGeo->endRate[i] == 0) {
581  myGeo->endRate[i] = myGeo->coarseRate[i];
582  ++myGeo->gCoarseNodesPerDir[i];
583  } else {
584  myGeo->gCoarseNodesPerDir[i] += 2;
585  }
586  } else {
587  myGeo->endRate[i] = 1;
588  myGeo->gCoarseNodesPerDir[i] = 1;
589  }
590  }
591 
592  myGeo->gNumCoarseNodes = myGeo->gCoarseNodesPerDir[0] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[2];
593 
594  for (LO i = 0; i < 3; ++i) {
595  if (i < myGeo->numDimensions) {
596  // Check whether the partition includes the "end" of the mesh which means that endRate will
597  // apply. Also make sure that endRate is not 0 which means that the mesh does not require a
598  // particular treatment at the boundaries.
599  if ((myGeo->startIndices[i] + myGeo->lFineNodesPerDir[i]) == myGeo->gFineNodesPerDir[i]) {
600  myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] - myGeo->endRate[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i] + 1;
601  if (myGeo->offsets[i] == 0) {
602  ++myGeo->lCoarseNodesPerDir[i];
603  }
604  } else {
605  myGeo->lCoarseNodesPerDir[i] = (myGeo->lFineNodesPerDir[i] + myGeo->offsets[i] - 1) / myGeo->coarseRate[i];
606  if (myGeo->offsets[i] == 0) {
607  ++myGeo->lCoarseNodesPerDir[i];
608  }
609  }
610  } else {
611  myGeo->lCoarseNodesPerDir[i] = 1;
612  }
613  // This would happen if the rank does not own any nodes but in that case a subcommunicator
614  // should be used so this should really not be a concern.
615  if (myGeo->lFineNodesPerDir[i] < 1) {
616  myGeo->lCoarseNodesPerDir[i] = 0;
617  }
618  }
619 
620  // Assuming linear interpolation, each fine point has contribution from 8 coarse points
621  // and each coarse point value gets injected.
622  // For systems of PDEs we assume that all dofs have the same P operator.
623  myGeo->lNumCoarseNodes = myGeo->lCoarseNodesPerDir[0] * myGeo->lCoarseNodesPerDir[1] * myGeo->lCoarseNodesPerDir[2];
624 
625  // For each direction, determine how many points (including ghosts) are required.
626  bool ghostedDir[6] = {false};
627  for (int dim = 0; dim < 3; ++dim) {
628  // The first branch of this if-statement will be used if the rank contains only one layer
629  // of nodes in direction i, that layer must also coincide with the boundary of the mesh
630  // and coarseRate[i] == endRate[i]...
631  if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1 &&
632  myGeo->startIndices[dim] % myGeo->coarseRate[dim] == 0) {
633  myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim] - 1;
634  } else {
635  myGeo->startGhostedCoarseNode[dim] = myGeo->startIndices[dim] / myGeo->coarseRate[dim];
636  }
637  myGeo->ghostedCoarseNodesPerDir[dim] = myGeo->lCoarseNodesPerDir[dim];
638  // Check whether face *low needs ghost nodes
639  if (myGeo->ghostInterface[2 * dim]) {
640  myGeo->ghostedCoarseNodesPerDir[dim] += 1;
641  ghostedDir[2 * dim] = true;
642  }
643  // Check whether face *hi needs ghost nodes
644  if (myGeo->ghostInterface[2 * dim + 1]) {
645  myGeo->ghostedCoarseNodesPerDir[dim] += 1;
646  ghostedDir[2 * dim + 1] = true;
647  }
648  }
649  myGeo->lNumGhostedNodes = myGeo->ghostedCoarseNodesPerDir[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0];
650  myGeo->lNumGhostNodes = myGeo->lNumGhostedNodes - myGeo->lNumCoarseNodes;
651  ghostedCoarseNodes->PIDs.resize(myGeo->lNumGhostedNodes);
652  ghostedCoarseNodes->LIDs.resize(myGeo->lNumGhostedNodes);
653  ghostedCoarseNodes->GIDs.resize(myGeo->lNumGhostedNodes);
654  ghostedCoarseNodes->coarseGIDs.resize(myGeo->lNumGhostedNodes);
655  ghostedCoarseNodes->colInds.resize(myGeo->lNumGhostedNodes);
656  lCoarseNodesGIDs[0].resize(myGeo->lNumCoarseNodes);
657  lCoarseNodesGIDs[1].resize(myGeo->lNumCoarseNodes);
658 
659  // Now the tricky part starts, the coarse nodes / ghosted coarse nodes need to be imported.
660  // This requires finding what their GID on the fine mesh is. They need to be ordered
661  // lexicographically to allow for fast sweeps through the mesh.
662 
663  // We loop over all ghosted coarse nodes by increasing global lexicographic order
664  Array<LO> coarseNodeCoarseIndices(3), coarseNodeFineIndices(3), ijk(3);
665  LO currentIndex = -1, countCoarseNodes = 0;
666  for (ijk[2] = 0; ijk[2] < myGeo->ghostedCoarseNodesPerDir[2]; ++ijk[2]) {
667  for (ijk[1] = 0; ijk[1] < myGeo->ghostedCoarseNodesPerDir[1]; ++ijk[1]) {
668  for (ijk[0] = 0; ijk[0] < myGeo->ghostedCoarseNodesPerDir[0]; ++ijk[0]) {
669  currentIndex = ijk[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + ijk[1] * myGeo->ghostedCoarseNodesPerDir[0] + ijk[0];
670  coarseNodeCoarseIndices[0] = myGeo->startGhostedCoarseNode[0] + ijk[0];
671  coarseNodeFineIndices[0] = coarseNodeCoarseIndices[0] * myGeo->coarseRate[0];
672  if (coarseNodeFineIndices[0] > myGeo->gFineNodesPerDir[0] - 1) {
673  coarseNodeFineIndices[0] = myGeo->gFineNodesPerDir[0] - 1;
674  }
675  coarseNodeCoarseIndices[1] = myGeo->startGhostedCoarseNode[1] + ijk[1];
676  coarseNodeFineIndices[1] = coarseNodeCoarseIndices[1] * myGeo->coarseRate[1];
677  if (coarseNodeFineIndices[1] > myGeo->gFineNodesPerDir[1] - 1) {
678  coarseNodeFineIndices[1] = myGeo->gFineNodesPerDir[1] - 1;
679  }
680  coarseNodeCoarseIndices[2] = myGeo->startGhostedCoarseNode[2] + ijk[2];
681  coarseNodeFineIndices[2] = coarseNodeCoarseIndices[2] * myGeo->coarseRate[2];
682  if (coarseNodeFineIndices[2] > myGeo->gFineNodesPerDir[2] - 1) {
683  coarseNodeFineIndices[2] = myGeo->gFineNodesPerDir[2] - 1;
684  }
685  GO myGID = 0, myCoarseGID = -1;
686  GO factor[3] = {};
687  factor[2] = myGeo->gNumFineNodes10;
688  factor[1] = myGeo->gFineNodesPerDir[0];
689  factor[0] = 1;
690  for (int dim = 0; dim < 3; ++dim) {
691  if (dim < myGeo->numDimensions) {
692  if (myGeo->startIndices[dim] - myGeo->offsets[dim] + ijk[dim] * myGeo->coarseRate[dim] < myGeo->gFineNodesPerDir[dim] - 1) {
693  myGID += (myGeo->startIndices[dim] - myGeo->offsets[dim] + ijk[dim] * myGeo->coarseRate[dim]) * factor[dim];
694  } else {
695  myGID += (myGeo->startIndices[dim] - myGeo->offsets[dim] + (ijk[dim] - 1) * myGeo->coarseRate[dim] + myGeo->endRate[dim]) * factor[dim];
696  }
697  }
698  }
699  myCoarseGID = coarseNodeCoarseIndices[0] + coarseNodeCoarseIndices[1] * myGeo->gCoarseNodesPerDir[0] + coarseNodeCoarseIndices[2] * myGeo->gCoarseNodesPerDir[1] * myGeo->gCoarseNodesPerDir[0];
700  ghostedCoarseNodes->GIDs[currentIndex] = myGID;
701  ghostedCoarseNodes->coarseGIDs[currentIndex] = myCoarseGID;
702  if ((!ghostedDir[0] || ijk[0] != 0) && (!ghostedDir[2] || ijk[1] != 0) && (!ghostedDir[4] || ijk[2] != 0) && (!ghostedDir[1] || ijk[0] != myGeo->ghostedCoarseNodesPerDir[0] - 1) && (!ghostedDir[3] || ijk[1] != myGeo->ghostedCoarseNodesPerDir[1] - 1) && (!ghostedDir[5] || ijk[2] != myGeo->ghostedCoarseNodesPerDir[2] - 1)) {
703  lCoarseNodesGIDs[0][countCoarseNodes] = myCoarseGID;
704  lCoarseNodesGIDs[1][countCoarseNodes] = myGID;
705  ++countCoarseNodes;
706  }
707  }
708  }
709  }
710  Array<int> ghostsPIDs(myGeo->lNumGhostedNodes);
711  Array<LO> ghostsLIDs(myGeo->lNumGhostedNodes);
712  fineCoordsMap->getRemoteIndexList(ghostedCoarseNodes->GIDs(),
713  ghostedCoarseNodes->PIDs(),
714  ghostedCoarseNodes->LIDs());
715 } // End GetCoarsePoint
716 
717 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
718 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
719  MakeGeneralGeometricP(RCP<GeometricData> myGeo,
720  const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, Node> >& fineCoords,
721  const LO nnzP, const LO dofsPerNode,
722  RCP<const Map>& stridedDomainMapP, RCP<Matrix>& Amat, RCP<Matrix>& P,
723  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, Node> >& coarseCoords,
724  RCP<NodesIDs> ghostedCoarseNodes, Array<Array<GO> > coarseNodesGIDs,
725  int interpolationOrder) const {
726  /* On termination, return the number of local prolongator columns owned by
727  * this processor.
728  *
729  * Input
730  * =====
731  * nNodes Number of fine level Blk Rows owned by this processor
732  * coarseRate Rate of coarsening in each spatial direction.
733  * endRate Rate of coarsening in each spatial direction for the last
734  * nodes in the mesh where an adaptive coarsening rate is
735  * required.
736  * nTerms Number of nonzero entries in the prolongation matrix.
737  * dofsPerNode Number of degrees-of-freedom per mesh node.
738  *
739  * Output
740  * =====
741  * So far nothing...
742  */
743 
744  using xdMV = Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>;
745  Xpetra::global_size_t OTI = Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid();
746 
747  LO myRank = Amat->getRowMap()->getComm()->getRank();
748  GO numGloCols = dofsPerNode * myGeo->gNumCoarseNodes;
749 
750  // Build maps necessary to create and fill complete the prolongator
751  // note: rowMapP == rangeMapP and colMapP != domainMapP
752  RCP<const Map> rowMapP = Amat->getDomainMap();
753 
754  RCP<const Map> domainMapP;
755 
756  RCP<const Map> colMapP;
757 
758  // At this point we need to create the column map which is a delicate operation.
759  // The entries in that map need to be ordered as follows:
760  // 1) first owned entries ordered by LID
761  // 2) second order the remaining entries by PID
762  // 3) entries with the same remote PID are ordered by GID.
763  // One piece of good news: myGeo->lNumCoarseNodes is the number of ownedNodes and
764  // myGeo->lNumGhostNodes is the number of remote nodes. The sorting can be limited to remote
765  // nodes as the owned ones are alreaded ordered by LID!
766 
767  {
768  std::vector<NodeID> colMapOrdering(myGeo->lNumGhostedNodes);
769  for (LO ind = 0; ind < myGeo->lNumGhostedNodes; ++ind) {
770  colMapOrdering[ind].GID = ghostedCoarseNodes->GIDs[ind];
771  if (ghostedCoarseNodes->PIDs[ind] == myRank) {
772  colMapOrdering[ind].PID = -1;
773  } else {
774  colMapOrdering[ind].PID = ghostedCoarseNodes->PIDs[ind];
775  }
776  colMapOrdering[ind].LID = ghostedCoarseNodes->LIDs[ind];
777  colMapOrdering[ind].lexiInd = ind;
778  }
779  std::sort(colMapOrdering.begin(), colMapOrdering.end(),
780  [](NodeID a, NodeID b) -> bool {
781  return ((a.PID < b.PID) || ((a.PID == b.PID) && (a.GID < b.GID)));
782  });
783 
784  Array<GO> colGIDs(dofsPerNode * myGeo->lNumGhostedNodes);
785  for (LO ind = 0; ind < myGeo->lNumGhostedNodes; ++ind) {
786  // Save the permutation calculated to go from Lexicographic indexing to column map indexing
787  ghostedCoarseNodes->colInds[colMapOrdering[ind].lexiInd] = ind;
788  for (LO dof = 0; dof < dofsPerNode; ++dof) {
789  colGIDs[dofsPerNode * ind + dof] = dofsPerNode * colMapOrdering[ind].GID + dof;
790  }
791  }
792  domainMapP = Xpetra::MapFactory<LO, GO, NO>::Build(rowMapP->lib(),
793  numGloCols,
794  colGIDs.view(0, dofsPerNode *
795  myGeo->lNumCoarseNodes),
796  rowMapP->getIndexBase(),
797  rowMapP->getComm());
798  colMapP = Xpetra::MapFactory<LO, GO, NO>::Build(rowMapP->lib(),
799  OTI,
800  colGIDs.view(0, colGIDs.size()),
801  rowMapP->getIndexBase(),
802  rowMapP->getComm());
803  } // End of scope for colMapOrdering and colGIDs
804 
805  std::vector<size_t> strideInfo(1);
806  strideInfo[0] = dofsPerNode;
807  stridedDomainMapP = Xpetra::StridedMapFactory<LO, GO, NO>::Build(domainMapP, strideInfo);
808 
809  // Build the map for the coarse level coordinates, create the associated MultiVector and fill it
810  // with an import from the fine coordinates MultiVector. As data is local this should not create
811  // communications during the importer creation.
812  RCP<const Map> coarseCoordsMap = MapFactory::Build(fineCoords->getMap()->lib(),
813  myGeo->gNumCoarseNodes,
814  coarseNodesGIDs[0](),
815  fineCoords->getMap()->getIndexBase(),
816  fineCoords->getMap()->getComm());
817  RCP<const Map> coarseCoordsFineMap = MapFactory::Build(fineCoords->getMap()->lib(),
818  myGeo->gNumCoarseNodes,
819  coarseNodesGIDs[1](),
820  fineCoords->getMap()->getIndexBase(),
821  fineCoords->getMap()->getComm());
822 
823  RCP<const Import> coarseImporter = ImportFactory::Build(fineCoords->getMap(),
824  coarseCoordsFineMap);
825  coarseCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>::Build(coarseCoordsFineMap,
826  myGeo->numDimensions);
827  coarseCoords->doImport(*fineCoords, *coarseImporter, Xpetra::INSERT);
828  coarseCoords->replaceMap(coarseCoordsMap);
829 
830  // Do the actual import using the fineCoords->getMap()
831  RCP<const Map> ghostMap = Xpetra::MapFactory<LO, GO, NO>::Build(fineCoords->getMap()->lib(),
832  OTI,
833  ghostedCoarseNodes->GIDs(),
834  fineCoords->getMap()->getIndexBase(),
835  fineCoords->getMap()->getComm());
836  RCP<const Import> ghostImporter = ImportFactory::Build(fineCoords->getMap(), ghostMap);
837  RCP<xdMV> ghostCoords = Xpetra::MultiVectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>::Build(ghostMap,
838  myGeo->numDimensions);
839  ghostCoords->doImport(*fineCoords, *ghostImporter, Xpetra::INSERT);
840 
841  P = rcp(new CrsMatrixWrap(rowMapP, colMapP, 0));
842  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
843 
844  ArrayRCP<size_t> iaP;
845  ArrayRCP<LO> jaP;
846  ArrayRCP<SC> valP;
847 
848  PCrs->allocateAllValues(nnzP, iaP, jaP, valP);
849 
850  ArrayView<size_t> ia = iaP();
851  ArrayView<LO> ja = jaP();
852  ArrayView<SC> val = valP();
853  ia[0] = 0;
854 
855  Array<ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> > ghostedCoords(3);
856  {
857  ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> tmp(ghostCoords->getLocalLength(), 0.0);
858  for (int dim = 0; dim < 3; ++dim) {
859  if (dim < myGeo->numDimensions) {
860  ghostedCoords[dim] = ghostCoords->getDataNonConst(dim);
861  } else {
862  ghostedCoords[dim] = tmp;
863  }
864  }
865  }
866 
867  // Declaration and assignment of fineCoords which holds the coordinates of the fine nodes in 3D.
868  // To do so we pull the nD coordinates from fineCoords and pad the rest with zero vectors...
869  RCP<Xpetra::Vector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO> > zeros = Xpetra::VectorFactory<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LO, GO, NO>::Build(fineCoords->getMap(), true);
870  ArrayRCP<ArrayRCP<typename Teuchos::ScalarTraits<Scalar>::coordinateType> > lFineCoords(3);
871  for (int dim = 0; dim < 3; ++dim) {
872  if (dim < myGeo->numDimensions) {
873  lFineCoords[dim] = fineCoords->getDataNonConst(dim);
874  } else {
875  lFineCoords[dim] = zeros->getDataNonConst(0);
876  }
877  }
878 
879  GO tStencil = 0;
880  for (int currentIndex = 0; currentIndex < myGeo->lNumFineNodes; ++currentIndex) {
881  Array<GO> ghostedIndices(3), firstInterpolationIndices(3);
882  Array<LO> interpolationPIDs(8), interpolationLIDs(8), interpolationGIDs(8);
883  Array<Array<typename Teuchos::ScalarTraits<Scalar>::coordinateType> > interpolationCoords(9);
884  interpolationCoords[0].resize(3);
885  GO firstInterpolationNodeIndex;
886  int nStencil = 0;
887  for (int dim = 0; dim < 3; ++dim) {
888  interpolationCoords[0][dim] = lFineCoords[dim][currentIndex];
889  }
890 
891  // Compute the ghosted (i,j,k) of the current node, that assumes (I,J,K)=(0,0,0) to be
892  // associated with the first node in ghostCoords
893  { // Scope for tmp
894  ghostedIndices[2] = currentIndex / (myGeo->lFineNodesPerDir[1] * myGeo->lFineNodesPerDir[0]);
895  LO tmp = currentIndex % (myGeo->lFineNodesPerDir[1] * myGeo->lFineNodesPerDir[0]);
896  ghostedIndices[1] = tmp / myGeo->lFineNodesPerDir[0];
897  ghostedIndices[0] = tmp % myGeo->lFineNodesPerDir[0];
898  for (int dim = 0; dim < 3; ++dim) {
899  ghostedIndices[dim] += myGeo->offsets[dim];
900  }
901  // A special case appears when the mesh is really coarse: it is possible for a rank to
902  // have a single coarse node in a given direction. If this happens on the highest i, j or k
903  // we need to "grab" a coarse node with a lower i, j, or k which leads us to add to the
904  // value of ghostedIndices
905  }
906  // No we find the indices of the coarse nodes used for interpolation simply by integer
907  // division.
908  for (int dim = 0; dim < 3; ++dim) {
909  firstInterpolationIndices[dim] = ghostedIndices[dim] / myGeo->coarseRate[dim];
910  // If you are on the edge of the local domain go back one coarse node, unless there is only
911  // one node on the local domain...
912  if (firstInterpolationIndices[dim] == myGeo->ghostedCoarseNodesPerDir[dim] - 1 && myGeo->ghostedCoarseNodesPerDir[dim] > 1) {
913  firstInterpolationIndices[dim] -= 1;
914  }
915  }
916  firstInterpolationNodeIndex =
917  firstInterpolationIndices[2] * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + firstInterpolationIndices[1] * myGeo->ghostedCoarseNodesPerDir[0] + firstInterpolationIndices[0];
918 
919  // We extract the coordinates and PIDs associated with each coarse node used during
920  // inteprolation in order to fill the prolongator correctly
921  {
922  LO ind = -1;
923  for (int k = 0; k < 2; ++k) {
924  for (int j = 0; j < 2; ++j) {
925  for (int i = 0; i < 2; ++i) {
926  ind = k * 4 + j * 2 + i;
927  interpolationLIDs[ind] = firstInterpolationNodeIndex + k * myGeo->ghostedCoarseNodesPerDir[1] * myGeo->ghostedCoarseNodesPerDir[0] + j * myGeo->ghostedCoarseNodesPerDir[0] + i;
928  if (ghostedCoarseNodes->PIDs[interpolationLIDs[ind]] == rowMapP->getComm()->getRank()) {
929  interpolationPIDs[ind] = -1;
930  } else {
931  interpolationPIDs[ind] = ghostedCoarseNodes->PIDs[interpolationLIDs[ind]];
932  }
933  interpolationGIDs[ind] = ghostedCoarseNodes->coarseGIDs[interpolationLIDs[ind]];
934 
935  interpolationCoords[ind + 1].resize(3);
936  for (int dim = 0; dim < 3; ++dim) {
937  interpolationCoords[ind + 1][dim] = ghostedCoords[dim][interpolationLIDs[ind]];
938  }
939  }
940  }
941  }
942  } // End of ind scope
943 
944  // Compute the actual geometric interpolation stencil
945  // LO stencilLength = static_cast<LO>(std::pow(interpolationOrder + 1, 3));
946  std::vector<double> stencil(8);
947  Array<GO> firstCoarseNodeFineIndices(3);
948  int rate[3] = {};
949  for (int dim = 0; dim < 3; ++dim) {
950  firstCoarseNodeFineIndices[dim] = firstInterpolationIndices[dim] * myGeo->coarseRate[dim];
951  if ((myGeo->startIndices[dim + 3] == myGeo->gFineNodesPerDir[dim] - 1) && (ghostedIndices[dim] >=
952  (myGeo->ghostedCoarseNodesPerDir[dim] - 2) * myGeo->coarseRate[dim])) {
953  rate[dim] = myGeo->endRate[dim];
954  } else {
955  rate[dim] = myGeo->coarseRate[dim];
956  }
957  }
958  Array<GO> trueGhostedIndices(3);
959  // This handles the case of a rank having a single node that also happens to be the last node
960  // in any direction. It might be more efficient to re-write the algo so that this is
961  // incorporated in the definition of ghostedIndices directly...
962  for (int dim = 0; dim < 3; ++dim) {
963  if (myGeo->startIndices[dim] == myGeo->gFineNodesPerDir[dim] - 1) {
964  trueGhostedIndices[dim] = ghostedIndices[dim] + rate[dim];
965  } else {
966  trueGhostedIndices[dim] = ghostedIndices[dim];
967  }
968  }
969  ComputeStencil(myGeo->numDimensions, trueGhostedIndices, firstCoarseNodeFineIndices, rate,
970  interpolationCoords, interpolationOrder, stencil);
971 
972  // Finally check whether the fine node is on a coarse: node, edge or face
973  // and select accordingly the non-zero values from the stencil and the
974  // corresponding column indices
975  Array<LO> nzIndStencil(8), permutation(8);
976  for (LO k = 0; k < 8; ++k) {
977  permutation[k] = k;
978  }
979  if (interpolationOrder == 0) {
980  nStencil = 1;
981  for (LO k = 0; k < 8; ++k) {
982  nzIndStencil[k] = static_cast<LO>(stencil[0]);
983  }
984  stencil[0] = 0.0;
985  stencil[nzIndStencil[0]] = 1.0;
986  } else if (interpolationOrder == 1) {
987  Array<GO> currentNodeGlobalFineIndices(3);
988  for (int dim = 0; dim < 3; ++dim) {
989  currentNodeGlobalFineIndices[dim] = ghostedIndices[dim] - myGeo->offsets[dim] + myGeo->startIndices[dim];
990  }
991  if (((ghostedIndices[0] % myGeo->coarseRate[0] == 0) || currentNodeGlobalFineIndices[0] == myGeo->gFineNodesPerDir[0] - 1) && ((ghostedIndices[1] % myGeo->coarseRate[1] == 0) || currentNodeGlobalFineIndices[1] == myGeo->gFineNodesPerDir[1] - 1) && ((ghostedIndices[2] % myGeo->coarseRate[2] == 0) || currentNodeGlobalFineIndices[2] == myGeo->gFineNodesPerDir[2] - 1)) {
992  if ((currentNodeGlobalFineIndices[0] == myGeo->gFineNodesPerDir[0] - 1) ||
993  (ghostedIndices[0] / myGeo->coarseRate[0] == myGeo->ghostedCoarseNodesPerDir[0] - 1)) {
994  nzIndStencil[0] += 1;
995  }
996  if (((currentNodeGlobalFineIndices[1] == myGeo->gFineNodesPerDir[1] - 1) ||
997  (ghostedIndices[1] / myGeo->coarseRate[1] == myGeo->ghostedCoarseNodesPerDir[1] - 1)) &&
998  (myGeo->numDimensions > 1)) {
999  nzIndStencil[0] += 2;
1000  }
1001  if (((currentNodeGlobalFineIndices[2] == myGeo->gFineNodesPerDir[2] - 1) ||
1002  (ghostedIndices[2] / myGeo->coarseRate[2] == myGeo->ghostedCoarseNodesPerDir[2] - 1)) &&
1003  myGeo->numDimensions > 2) {
1004  nzIndStencil[0] += 4;
1005  }
1006  nStencil = 1;
1007  for (LO k = 0; k < 8; ++k) {
1008  nzIndStencil[k] = nzIndStencil[0];
1009  }
1010  } else {
1011  nStencil = 8;
1012  for (LO k = 0; k < 8; ++k) {
1013  nzIndStencil[k] = k;
1014  }
1015  }
1016  }
1017 
1018  // Here the values are filled in the Crs matrix arrays
1019  // This is basically the only place these variables are modified
1020  // Hopefully this makes handling system of PDEs easy!
1021 
1022  // Loop on dofsPerNode and process each row for the current Node
1023 
1024  // Sort nodes by PIDs using stable sort to keep sublist ordered by LIDs and GIDs
1025  sh_sort2(interpolationPIDs.begin(), interpolationPIDs.end(),
1026  permutation.begin(), permutation.end());
1027 
1028  GO colInd;
1029  for (LO j = 0; j < dofsPerNode; ++j) {
1030  ia[currentIndex * dofsPerNode + j + 1] = ia[currentIndex * dofsPerNode + j] + nStencil;
1031  for (LO k = 0; k < nStencil; ++k) {
1032  colInd = ghostedCoarseNodes->colInds[interpolationLIDs[nzIndStencil[permutation[k]]]];
1033  ja[ia[currentIndex * dofsPerNode + j] + k] = colInd * dofsPerNode + j;
1034  val[ia[currentIndex * dofsPerNode + j] + k] = stencil[nzIndStencil[permutation[k]]];
1035  }
1036  // Add the stencil for each degree of freedom.
1037  tStencil += nStencil;
1038  }
1039  } // End loop over fine nodes
1040 
1041  if (rowMapP->lib() == Xpetra::UseTpetra) {
1042  // - Cannot resize for Epetra, as it checks for same pointers
1043  // - Need to resize for Tpetra, as it check ().size() == ia[numRows]
1044  // NOTE: these invalidate ja and val views
1045  jaP.resize(tStencil);
1046  valP.resize(tStencil);
1047  }
1048 
1049  // Set the values of the prolongation operators into the CrsMatrix P and call FillComplete
1050  PCrs->setAllValues(iaP, jaP, valP);
1051  PCrs->expertStaticFillComplete(domainMapP, rowMapP);
1052 } // End MakeGeneralGeometricP
1053 
1054 // template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1055 // void BlackBoxPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetGeometricData(
1056 // RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType,LO,GO,NO> >& coordinates, const Array<LO> coarseRate,
1057 // const Array<GO> gFineNodesPerDir, const Array<LO> lFineNodesPerDir, const LO BlkSize,
1058 // Array<GO>& gIndices, Array<LO>& myOffset, Array<bool>& ghostInterface, Array<LO>& endRate,
1059 // Array<GO>& gCoarseNodesPerDir, Array<LO>& lCoarseNodesPerDir, Array<GO>& ghostGIDs,
1060 // Array<GO>& coarseNodesGIDs, Array<GO>& colGIDs, GO& gNumCoarseNodes, LO& lNumCoarseNodes,
1061 // ArrayRCP<Array<double> > coarseNodes) const {
1062 
1063 // RCP<const Map> coordinatesMap = coordinates->getMap();
1064 // LO numDimensions = coordinates->getNumVectors();
1065 
1066 // // Using the coarsening rate and the fine level data,
1067 // // compute coarse level data
1068 
1069 // // Phase 1 //
1070 // // ------------------------------------------------------------------ //
1071 // // We first start by finding small informations on the mesh such as //
1072 // // the number of coarse nodes (local and global) and the number of //
1073 // // ghost nodes / the end rate of coarsening. //
1074 // // ------------------------------------------------------------------ //
1075 // GO minGlobalIndex = coordinatesMap->getMinGlobalIndex();
1076 // {
1077 // GO tmp;
1078 // gIndices[2] = minGlobalIndex / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
1079 // tmp = minGlobalIndex % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
1080 // gIndices[1] = tmp / gFineNodesPerDir[0];
1081 // gIndices[0] = tmp % gFineNodesPerDir[0];
1082 
1083 // myOffset[2] = gIndices[2] % coarseRate[2];
1084 // myOffset[1] = gIndices[1] % coarseRate[1];
1085 // myOffset[0] = gIndices[0] % coarseRate[0];
1086 // }
1087 
1088 // // Check whether ghost nodes are needed in each direction
1089 // for(LO i=0; i < numDimensions; ++i) {
1090 // if((gIndices[i] != 0) && (gIndices[i] % coarseRate[i] > 0)) {
1091 // ghostInterface[2*i] = true;
1092 // }
1093 // if(((gIndices[i] + lFineNodesPerDir[i]) != gFineNodesPerDir[i]) && ((gIndices[i] + lFineNodesPerDir[i] - 1) % coarseRate[i] > 0)) {
1094 // ghostInterface[2*i + 1] = true;
1095 // }
1096 // }
1097 
1098 // for(LO i = 0; i < 3; ++i) {
1099 // if(i < numDimensions) {
1100 // lCoarseNodesPerDir[i] = (lFineNodesPerDir[i] + myOffset[i] - 1) / coarseRate[i];
1101 // if(myOffset[i] == 0) { ++lCoarseNodesPerDir[i]; }
1102 // gCoarseNodesPerDir[i] = (gFineNodesPerDir[i] - 1) / coarseRate[i];
1103 // endRate[i] = (gFineNodesPerDir[i] - 1) % coarseRate[i];
1104 // if(endRate[i] == 0) {
1105 // ++gCoarseNodesPerDir[i];
1106 // endRate[i] = coarseRate[i];
1107 // }
1108 // } else {
1109 // // Most quantities need to be set to 1 for extra dimensions
1110 // // this is rather logical, an x-y plane is like a single layer
1111 // // of nodes in the z direction...
1112 // gCoarseNodesPerDir[i] = 1;
1113 // lCoarseNodesPerDir[i] = 1;
1114 // endRate[i] = 1;
1115 // }
1116 // }
1117 
1118 // gNumCoarseNodes = gCoarseNodesPerDir[0]*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[2];
1119 // lNumCoarseNodes = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];
1120 
1121 // // For each direction, determine how many ghost points are required.
1122 // LO lNumGhostNodes = 0;
1123 // {
1124 // const int complementaryIndices[3][2] = {{1,2}, {0,2}, {0,1}};
1125 // LO tmp = 0;
1126 // for(LO i = 0; i < 3; ++i) {
1127 // // Check whether a face in direction i needs ghost nodes
1128 // if(ghostInterface[2*i] || ghostInterface[2*i+1]) {
1129 // if(i == 0) {tmp = lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];}
1130 // if(i == 1) {tmp = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[2];}
1131 // if(i == 2) {tmp = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1];}
1132 // }
1133 // // If both faces in direction i need nodes, double the number of ghost nodes
1134 // if(ghostInterface[2*i] && ghostInterface[2*i+1]) {tmp = 2*tmp;}
1135 // lNumGhostNodes += tmp;
1136 
1137 // // The corners and edges need to be checked in 2D / 3D to add more ghosts nodes
1138 // for(LO j = 0; j < 2; ++j) {
1139 // for(LO k = 0; k < 2; ++k) {
1140 // // Check if two adjoining faces need ghost nodes and then add their common edge
1141 // if(ghostInterface[2*complementaryIndices[i][0]+j] && ghostInterface[2*complementaryIndices[i][1]+k]) {
1142 // lNumGhostNodes += lCoarseNodesPerDir[i];
1143 // // Add corners if three adjoining faces need ghost nodes,
1144 // // but add them only once! Hence when i == 0.
1145 // if(ghostInterface[2*i] && (i == 0)) { lNumGhostNodes += 1; }
1146 // if(ghostInterface[2*i+1] && (i == 0)) { lNumGhostNodes += 1; }
1147 // }
1148 // }
1149 // }
1150 // tmp = 0;
1151 // }
1152 // } // end of scope for tmp and complementaryIndices
1153 
1154 // // Phase 2 //
1155 // // ------------------------------------------------------------------ //
1156 // // Build a list of GIDs to import the required ghost nodes. //
1157 // // The ordering of the ghosts nodes will be as natural as possible, //
1158 // // i.e. it should follow the GID ordering of the mesh. //
1159 // // ------------------------------------------------------------------ //
1160 
1161 // // Saddly we have to more or less redo what was just done to figure out the value of lNumGhostNodes,
1162 // // there might be some optimization possibility here...
1163 // ghostGIDs.resize(lNumGhostNodes);
1164 // LO countGhosts = 0;
1165 // // Get the GID of the first point on the current partition.
1166 // GO startingGID = minGlobalIndex;
1167 // Array<GO> startingIndices(3);
1168 // // We still want ghost nodes even if have with a 0 offset,
1169 // // except when we are on a boundary
1170 // if(ghostInterface[4] && (myOffset[2] == 0)) {
1171 // if(gIndices[2] + coarseRate[2] > gFineNodesPerDir[2]) {
1172 // startingGID -= endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1173 // } else {
1174 // startingGID -= coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1175 // }
1176 // } else {
1177 // startingGID -= myOffset[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1178 // }
1179 // if(ghostInterface[2] && (myOffset[1] == 0)) {
1180 // if(gIndices[1] + coarseRate[1] > gFineNodesPerDir[1]) {
1181 // startingGID -= endRate[1]*gFineNodesPerDir[0];
1182 // } else {
1183 // startingGID -= coarseRate[1]*gFineNodesPerDir[0];
1184 // }
1185 // } else {
1186 // startingGID -= myOffset[1]*gFineNodesPerDir[0];
1187 // }
1188 // if(ghostInterface[0] && (myOffset[0] == 0)) {
1189 // if(gIndices[0] + coarseRate[0] > gFineNodesPerDir[0]) {
1190 // startingGID -= endRate[0];
1191 // } else {
1192 // startingGID -= coarseRate[0];
1193 // }
1194 // } else {
1195 // startingGID -= myOffset[0];
1196 // }
1197 
1198 // { // scope for tmp
1199 // GO tmp;
1200 // startingIndices[2] = startingGID / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
1201 // tmp = startingGID % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
1202 // startingIndices[1] = tmp / gFineNodesPerDir[0];
1203 // startingIndices[0] = tmp % gFineNodesPerDir[0];
1204 // }
1205 
1206 // GO ghostOffset[3] = {0, 0, 0};
1207 // LO lengthZero = lCoarseNodesPerDir[0];
1208 // LO lengthOne = lCoarseNodesPerDir[1];
1209 // LO lengthTwo = lCoarseNodesPerDir[2];
1210 // if(ghostInterface[0]) {++lengthZero;}
1211 // if(ghostInterface[1]) {++lengthZero;}
1212 // if(ghostInterface[2]) {++lengthOne;}
1213 // if(ghostInterface[3]) {++lengthOne;}
1214 // if(ghostInterface[4]) {++lengthTwo;}
1215 // if(ghostInterface[5]) {++lengthTwo;}
1216 
1217 // // First check the bottom face as it will have the lowest GIDs
1218 // if(ghostInterface[4]) {
1219 // ghostOffset[2] = startingGID;
1220 // for(LO j = 0; j < lengthOne; ++j) {
1221 // if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1222 // ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
1223 // } else {
1224 // ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1225 // }
1226 // for(LO k = 0; k < lengthZero; ++k) {
1227 // if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1228 // ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1229 // } else {
1230 // ghostOffset[0] = k*coarseRate[0];
1231 // }
1232 // // If the partition includes a changed rate at the edge, ghost nodes need to be picked carefully.
1233 // // This if statement is repeated each time a ghost node is selected.
1234 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1235 // ++countGhosts;
1236 // }
1237 // }
1238 // }
1239 
1240 // // Sweep over the lCoarseNodesPerDir[2] coarse layers in direction 2 and gather necessary ghost nodes
1241 // // located on these layers.
1242 // for(LO i = 0; i < lengthTwo; ++i) {
1243 // // Exclude the cases where ghost nodes exists on the faces in directions 2, these faces are swept
1244 // // seperatly for efficiency.
1245 // if( !((i == lengthTwo-1) && ghostInterface[5]) && !((i == 0) && ghostInterface[4]) ) {
1246 // // Set the ghostOffset in direction 2 taking into account a possible endRate different
1247 // // from the regular coarseRate.
1248 // if( (i == lengthTwo-1) && !ghostInterface[5] ) {
1249 // ghostOffset[2] = startingGID + ((i-1)*coarseRate[2] + endRate[2])*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1250 // } else {
1251 // ghostOffset[2] = startingGID + i*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1252 // }
1253 // for(LO j = 0; j < lengthOne; ++j) {
1254 // if( (j == 0) && ghostInterface[2] ) {
1255 // for(LO k = 0; k < lengthZero; ++k) {
1256 // if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1257 // if(k == 0) {
1258 // ghostOffset[0] = endRate[0];
1259 // } else {
1260 // ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1261 // }
1262 // } else {
1263 // ghostOffset[0] = k*coarseRate[0];
1264 // }
1265 // // In this case j == 0 so ghostOffset[1] is zero and can be ignored
1266 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[0];
1267 // ++countGhosts;
1268 // }
1269 // } else if( (j == lengthOne-1) && ghostInterface[3] ) {
1270 // // Set the ghostOffset in direction 1 taking into account a possible endRate different
1271 // // from the regular coarseRate.
1272 // if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1273 // ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
1274 // } else {
1275 // ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1276 // }
1277 // for(LO k = 0; k < lengthZero; ++k) {
1278 // if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1279 // ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1280 // } else {
1281 // ghostOffset[0] = k*coarseRate[0];
1282 // }
1283 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1284 // ++countGhosts;
1285 // }
1286 // } else {
1287 // // Set ghostOffset[1] for side faces sweep
1288 // if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
1289 // ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
1290 // } else {
1291 // ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1292 // }
1293 
1294 // // Set ghostOffset[0], ghostGIDs and countGhosts
1295 // if(ghostInterface[0]) { // In that case ghostOffset[0]==0, so we can ignore it
1296 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1];
1297 // ++countGhosts;
1298 // }
1299 // if(ghostInterface[1]) { // Grab ghost point at the end of direction 0.
1300 // if( (startingIndices[0] + (lengthZero-1)*coarseRate[0]) > gFineNodesPerDir[0] - 1 ) {
1301 // if(lengthZero > 2) {
1302 // ghostOffset[0] = (lengthZero-2)*coarseRate[0] + endRate[0];
1303 // } else {
1304 // ghostOffset[0] = endRate[0];
1305 // }
1306 // } else {
1307 // ghostOffset[0] = (lengthZero-1)*coarseRate[0];
1308 // }
1309 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1310 // ++countGhosts;
1311 // }
1312 // }
1313 // }
1314 // }
1315 // }
1316 
1317 // // Finally check the top face
1318 // if(ghostInterface[5]) {
1319 // if( startingIndices[2] + (lengthTwo-1)*coarseRate[2] + 1 > gFineNodesPerDir[2] ) {
1320 // ghostOffset[2] = startingGID + ((lengthTwo-2)*coarseRate[2] + endRate[2])*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1321 // } else {
1322 // ghostOffset[2] = startingGID + (lengthTwo-1)*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1323 // }
1324 // for(LO j = 0; j < lengthOne; ++j) {
1325 // if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) { // && !ghostInterface[3] ) {
1326 // ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
1327 // } else {
1328 // ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
1329 // }
1330 // for(LO k = 0; k < lengthZero; ++k) {
1331 // if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
1332 // ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
1333 // } else {
1334 // ghostOffset[0] = k*coarseRate[0];
1335 // }
1336 // ghostGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
1337 // ++countGhosts;
1338 // }
1339 // }
1340 // }
1341 
1342 // // Phase 3 //
1343 // // ------------------------------------------------------------------ //
1344 // // Final phase of this function, lists are being built to create the //
1345 // // column and domain maps of the projection as well as the map and //
1346 // // arrays for the coarse coordinates multivector. //
1347 // // ------------------------------------------------------------------ //
1348 
1349 // Array<GO> gCoarseNodesGIDs(lNumCoarseNodes);
1350 // LO currentNode, offset2, offset1, offset0;
1351 // // Find the GIDs of the coarse nodes on the partition.
1352 // for(LO ind2 = 0; ind2 < lCoarseNodesPerDir[2]; ++ind2) {
1353 // if(myOffset[2] == 0) {
1354 // offset2 = startingIndices[2] + myOffset[2];
1355 // } else {
1356 // if(startingIndices[2] + endRate[2] == gFineNodesPerDir[2] - 1) {
1357 // offset2 = startingIndices[2] + endRate[2];
1358 // } else {
1359 // offset2 = startingIndices[2] + coarseRate[2];
1360 // }
1361 // }
1362 // if(offset2 + ind2*coarseRate[2] > gFineNodesPerDir[2] - 1) {
1363 // offset2 += (ind2 - 1)*coarseRate[2] + endRate[2];
1364 // } else {
1365 // offset2 += ind2*coarseRate[2];
1366 // }
1367 // offset2 = offset2*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1368 
1369 // for(LO ind1 = 0; ind1 < lCoarseNodesPerDir[1]; ++ind1) {
1370 // if(myOffset[1] == 0) {
1371 // offset1 = startingIndices[1] + myOffset[1];
1372 // } else {
1373 // if(startingIndices[1] + endRate[1] == gFineNodesPerDir[1] - 1) {
1374 // offset1 = startingIndices[1] + endRate[1];
1375 // } else {
1376 // offset1 = startingIndices[1] + coarseRate[1];
1377 // }
1378 // }
1379 // if(offset1 + ind1*coarseRate[1] > gFineNodesPerDir[1] - 1) {
1380 // offset1 += (ind1 - 1)*coarseRate[1] + endRate[1];
1381 // } else {
1382 // offset1 += ind1*coarseRate[1];
1383 // }
1384 // offset1 = offset1*gFineNodesPerDir[0];
1385 // for(LO ind0 = 0; ind0 < lCoarseNodesPerDir[0]; ++ind0) {
1386 // offset0 = startingIndices[0];
1387 // if(myOffset[0] == 0) {
1388 // offset0 += ind0*coarseRate[0];
1389 // } else {
1390 // offset0 += (ind0 + 1)*coarseRate[0];
1391 // }
1392 // if(offset0 > gFineNodesPerDir[0] - 1) {offset0 += endRate[0] - coarseRate[0];}
1393 
1394 // currentNode = ind2*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]
1395 // + ind1*lCoarseNodesPerDir[0]
1396 // + ind0;
1397 // gCoarseNodesGIDs[currentNode] = offset2 + offset1 + offset0;
1398 // }
1399 // }
1400 // }
1401 
1402 // // Actual loop over all the coarse/ghost nodes to find their index on the coarse mesh
1403 // // and the corresponding dofs that will need to be added to colMapP.
1404 // colGIDs.resize(BlkSize*(lNumCoarseNodes+lNumGhostNodes));
1405 // coarseNodesGIDs.resize(lNumCoarseNodes);
1406 // for(LO i = 0; i < numDimensions; ++i) {coarseNodes[i].resize(lNumCoarseNodes);}
1407 // GO fineNodesPerCoarseSlab = coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1408 // GO fineNodesEndCoarseSlab = endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
1409 // GO fineNodesPerCoarsePlane = coarseRate[1]*gFineNodesPerDir[0];
1410 // GO fineNodesEndCoarsePlane = endRate[1]*gFineNodesPerDir[0];
1411 // GO coarseNodesPerCoarseLayer = gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0];
1412 // GO gCoarseNodeOnCoarseGridGID;
1413 // LO gInd[3], lCol;
1414 // Array<int> ghostPIDs (lNumGhostNodes);
1415 // Array<LO> ghostLIDs (lNumGhostNodes);
1416 // Array<LO> ghostPermut(lNumGhostNodes);
1417 // for(LO k = 0; k < lNumGhostNodes; ++k) {ghostPermut[k] = k;}
1418 // coordinatesMap->getRemoteIndexList(ghostGIDs, ghostPIDs, ghostLIDs);
1419 // sh_sort_permute(ghostPIDs.begin(),ghostPIDs.end(), ghostPermut.begin(),ghostPermut.end());
1420 
1421 // { // scope for tmpInds, tmpVars and tmp.
1422 // GO tmpInds[3], tmpVars[2];
1423 // LO tmp;
1424 // // Loop over the coarse nodes of the partition and add them to colGIDs
1425 // // that will be used to construct the column and domain maps of P as well
1426 // // as to construct the coarse coordinates map.
1427 // // for(LO col = 0; col < lNumCoarseNodes; ++col) { // This should most likely be replaced by loops of lCoarseNodesPerDir[] to simplify arithmetics
1428 // LO col = 0;
1429 // LO firstCoarseNodeInds[3], currentCoarseNode;
1430 // for(LO dim = 0; dim < 3; ++dim) {
1431 // if(myOffset[dim] == 0) {
1432 // firstCoarseNodeInds[dim] = 0;
1433 // } else {
1434 // firstCoarseNodeInds[dim] = coarseRate[dim] - myOffset[dim];
1435 // }
1436 // }
1437 // Array<ArrayRCP<const typename Teuchos::ScalarTraits<Scalar>::coordinateType> > fineNodes(numDimensions);
1438 // for(LO dim = 0; dim < numDimensions; ++dim) {fineNodes[dim] = coordinates->getData(dim);}
1439 // for(LO k = 0; k < lCoarseNodesPerDir[2]; ++k) {
1440 // for(LO j = 0; j < lCoarseNodesPerDir[1]; ++j) {
1441 // for(LO i = 0; i < lCoarseNodesPerDir[0]; ++i) {
1442 // col = k*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0] + j*lCoarseNodesPerDir[0] + i;
1443 
1444 // // Check for endRate
1445 // currentCoarseNode = 0;
1446 // if(firstCoarseNodeInds[0] + i*coarseRate[0] > lFineNodesPerDir[0] - 1) {
1447 // currentCoarseNode += firstCoarseNodeInds[0] + (i-1)*coarseRate[0] + endRate[0];
1448 // } else {
1449 // currentCoarseNode += firstCoarseNodeInds[0] + i*coarseRate[0];
1450 // }
1451 // if(firstCoarseNodeInds[1] + j*coarseRate[1] > lFineNodesPerDir[1] - 1) {
1452 // currentCoarseNode += (firstCoarseNodeInds[1] + (j-1)*coarseRate[1] + endRate[1])*lFineNodesPerDir[0];
1453 // } else {
1454 // currentCoarseNode += (firstCoarseNodeInds[1] + j*coarseRate[1])*lFineNodesPerDir[0];
1455 // }
1456 // if(firstCoarseNodeInds[2] + k*coarseRate[2] > lFineNodesPerDir[2] - 1) {
1457 // currentCoarseNode += (firstCoarseNodeInds[2] + (k-1)*coarseRate[2] + endRate[2])*lFineNodesPerDir[1]*lFineNodesPerDir[0];
1458 // } else {
1459 // currentCoarseNode += (firstCoarseNodeInds[2] + k*coarseRate[2])*lFineNodesPerDir[1]*lFineNodesPerDir[0];
1460 // }
1461 // // Load coordinates
1462 // for(LO dim = 0; dim < numDimensions; ++dim) {
1463 // coarseNodes[dim][col] = fineNodes[dim][currentCoarseNode];
1464 // }
1465 
1466 // if((endRate[2] != coarseRate[2]) && (gCoarseNodesGIDs[col] > (gCoarseNodesPerDir[2] - 2)*fineNodesPerCoarseSlab + fineNodesEndCoarseSlab - 1)) {
1467 // tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab + 1;
1468 // tmpVars[0] = gCoarseNodesGIDs[col] - (tmpInds[2] - 1)*fineNodesPerCoarseSlab - fineNodesEndCoarseSlab;
1469 // } else {
1470 // tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab;
1471 // tmpVars[0] = gCoarseNodesGIDs[col] % fineNodesPerCoarseSlab;
1472 // }
1473 // if((endRate[1] != coarseRate[1]) && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane + fineNodesEndCoarsePlane - 1)) {
1474 // tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
1475 // tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane - fineNodesEndCoarsePlane;
1476 // } else {
1477 // tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
1478 // tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
1479 // }
1480 // if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
1481 // tmpInds[0] = gCoarseNodesPerDir[0] - 1;
1482 // } else {
1483 // tmpInds[0] = tmpVars[1] / coarseRate[0];
1484 // }
1485 // gInd[2] = col / (lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]);
1486 // tmp = col % (lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]);
1487 // gInd[1] = tmp / lCoarseNodesPerDir[0];
1488 // gInd[0] = tmp % lCoarseNodesPerDir[0];
1489 // lCol = gInd[2]*(lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]) + gInd[1]*lCoarseNodesPerDir[0] + gInd[0];
1490 // gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
1491 // coarseNodesGIDs[lCol] = gCoarseNodeOnCoarseGridGID;
1492 // for(LO dof = 0; dof < BlkSize; ++dof) {
1493 // colGIDs[BlkSize*lCol + dof] = BlkSize*gCoarseNodeOnCoarseGridGID + dof;
1494 // }
1495 // }
1496 // }
1497 // }
1498 // // Now loop over the ghost nodes of the partition to add them to colGIDs
1499 // // since they will need to be included in the column map of P
1500 // for(col = lNumCoarseNodes; col < lNumCoarseNodes + lNumGhostNodes; ++col) {
1501 // if((endRate[2] != coarseRate[2]) && (ghostGIDs[ghostPermut[col - lNumCoarseNodes]] > (gCoarseNodesPerDir[2] - 2)*fineNodesPerCoarseSlab + fineNodesEndCoarseSlab - 1)) {
1502 // tmpInds[2] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab + 1;
1503 // tmpVars[0] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] - (tmpInds[2] - 1)*fineNodesPerCoarseSlab - fineNodesEndCoarseSlab;
1504 // } else {
1505 // tmpInds[2] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab;
1506 // tmpVars[0] = ghostGIDs[ghostPermut[col - lNumCoarseNodes]] % fineNodesPerCoarseSlab;
1507 // }
1508 // if((endRate[1] != coarseRate[1]) && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane + fineNodesEndCoarsePlane - 1)) {
1509 // tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
1510 // tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane - fineNodesEndCoarsePlane;
1511 // } else {
1512 // tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
1513 // tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
1514 // }
1515 // if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
1516 // tmpInds[0] = gCoarseNodesPerDir[0] - 1;
1517 // } else {
1518 // tmpInds[0] = tmpVars[1] / coarseRate[0];
1519 // }
1520 // gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
1521 // for(LO dof = 0; dof < BlkSize; ++dof) {
1522 // colGIDs[BlkSize*col + dof] = BlkSize*gCoarseNodeOnCoarseGridGID + dof;
1523 // }
1524 // }
1525 // } // End of scope for tmpInds, tmpVars and tmp
1526 
1527 // } // GetGeometricData()
1528 
1529 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1530 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeStencil(
1531  const LO numDimensions, const Array<GO> currentNodeIndices,
1532  const Array<GO> coarseNodeIndices, const LO rate[3],
1533  const Array<Array<typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coord, const int interpolationOrder,
1534  std::vector<double>& stencil) const {
1535  TEUCHOS_TEST_FOR_EXCEPTION((interpolationOrder > 1) || (interpolationOrder < 0),
1536  Exceptions::RuntimeError,
1537  "The interpolation order can be set to 0 or 1 only.");
1538 
1539  if (interpolationOrder == 0) {
1540  ComputeConstantInterpolationStencil(numDimensions, currentNodeIndices, coarseNodeIndices,
1541  rate, stencil);
1542  } else if (interpolationOrder == 1) {
1543  ComputeLinearInterpolationStencil(numDimensions, coord, stencil);
1544  }
1545 
1546 } // End ComputeStencil
1547 
1548 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1549 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1550  ComputeConstantInterpolationStencil(const LO numDimensions, const Array<GO> currentNodeIndices,
1551  const Array<GO> coarseNodeIndices, const LO rate[3],
1552  std::vector<double>& stencil) const {
1553  LO coarseNode = 0;
1554  if (numDimensions > 2) {
1555  if ((currentNodeIndices[2] - coarseNodeIndices[2]) > (rate[2] / 2)) {
1556  coarseNode += 4;
1557  }
1558  }
1559  if (numDimensions > 1) {
1560  if ((currentNodeIndices[1] - coarseNodeIndices[1]) > (rate[1] / 2)) {
1561  coarseNode += 2;
1562  }
1563  }
1564  if ((currentNodeIndices[0] - coarseNodeIndices[0]) > (rate[0] / 2)) {
1565  coarseNode += 1;
1566  }
1567  stencil[0] = coarseNode;
1568 
1569 } // ComputeConstantInterpolationStencil
1570 
1571 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1572 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1573  ComputeLinearInterpolationStencil(const LO numDimensions, const Array<Array<typename Teuchos::ScalarTraits<Scalar>::coordinateType> > coord,
1574  std::vector<double>& stencil)
1575  const {
1576  // 7 8 Find xi, eta and zeta such that
1577  // x---------x
1578  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
1579  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
1580  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
1581  // | | *P | |
1582  // | x------|--x We can do this with a Newton solver:
1583  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
1584  // |/ |/ We compute the Jacobian and iterate until convergence...
1585  // z y x---------x
1586  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
1587  // |/ give us the weights for the interpolation stencil!
1588  // o---x
1589  //
1590 
1591  Teuchos::SerialDenseMatrix<LO, double> Jacobian(numDimensions, numDimensions);
1592  Teuchos::SerialDenseVector<LO, double> residual(numDimensions);
1593  Teuchos::SerialDenseVector<LO, double> solutionDirection(numDimensions);
1594  Teuchos::SerialDenseVector<LO, double> paramCoords(numDimensions);
1595  Teuchos::SerialDenseSolver<LO, double> problem;
1596  LO numTerms = std::pow(2, numDimensions), iter = 0, max_iter = 5;
1597  double functions[4][8], norm_ref = 1, norm2 = 1, tol = 1e-5;
1598  paramCoords.size(numDimensions);
1599 
1600  while ((iter < max_iter) && (norm2 > tol * norm_ref)) {
1601  ++iter;
1602  norm2 = 0;
1603  solutionDirection.size(numDimensions);
1604  residual.size(numDimensions);
1605  Jacobian = 0.0;
1606 
1607  // Compute Jacobian and Residual
1608  GetInterpolationFunctions(numDimensions, paramCoords, functions);
1609  for (LO i = 0; i < numDimensions; ++i) {
1610  residual(i) = coord[0][i]; // Add coordinates from point of interest
1611  for (LO k = 0; k < numTerms; ++k) {
1612  residual(i) -= functions[0][k] * coord[k + 1][i]; // Remove contribution from all coarse points
1613  }
1614  if (iter == 1) {
1615  norm_ref += residual(i) * residual(i);
1616  if (i == numDimensions - 1) {
1617  norm_ref = std::sqrt(norm_ref);
1618  }
1619  }
1620 
1621  for (LO j = 0; j < numDimensions; ++j) {
1622  for (LO k = 0; k < numTerms; ++k) {
1623  Jacobian(i, j) += functions[j + 1][k] * coord[k + 1][i];
1624  }
1625  }
1626  }
1627 
1628  // Set Jacobian, Vectors and solve problem
1629  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
1630  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
1631  problem.factorWithEquilibration(true);
1632  problem.solve();
1633  problem.unequilibrateLHS();
1634 
1635  for (LO i = 0; i < numDimensions; ++i) {
1636  paramCoords(i) = paramCoords(i) + solutionDirection(i);
1637  }
1638 
1639  // Recompute Residual norm
1640  GetInterpolationFunctions(numDimensions, paramCoords, functions);
1641  for (LO i = 0; i < numDimensions; ++i) {
1642  double tmp = coord[0][i];
1643  for (LO k = 0; k < numTerms; ++k) {
1644  tmp -= functions[0][k] * coord[k + 1][i];
1645  }
1646  norm2 += tmp * tmp;
1647  tmp = 0;
1648  }
1649  norm2 = std::sqrt(norm2);
1650  }
1651 
1652  // Load the interpolation values onto the stencil.
1653  for (LO i = 0; i < 8; ++i) {
1654  stencil[i] = functions[0][i];
1655  }
1656 
1657 } // End ComputeLinearInterpolationStencil
1658 
1659 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1660 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1661  GetInterpolationFunctions(const LO numDimensions,
1662  const Teuchos::SerialDenseVector<LO, double> parameters,
1663  double functions[4][8]) const {
1664  double xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
1665  if (numDimensions == 1) {
1666  xi = parameters[0];
1667  denominator = 2.0;
1668  } else if (numDimensions == 2) {
1669  xi = parameters[0];
1670  eta = parameters[1];
1671  denominator = 4.0;
1672  } else if (numDimensions == 3) {
1673  xi = parameters[0];
1674  eta = parameters[1];
1675  zeta = parameters[2];
1676  denominator = 8.0;
1677  }
1678 
1679  functions[0][0] = (1.0 - xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
1680  functions[0][1] = (1.0 + xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
1681  functions[0][2] = (1.0 - xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
1682  functions[0][3] = (1.0 + xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
1683  functions[0][4] = (1.0 - xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
1684  functions[0][5] = (1.0 + xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
1685  functions[0][6] = (1.0 - xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
1686  functions[0][7] = (1.0 + xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
1687 
1688  functions[1][0] = -(1.0 - eta) * (1.0 - zeta) / denominator;
1689  functions[1][1] = (1.0 - eta) * (1.0 - zeta) / denominator;
1690  functions[1][2] = -(1.0 + eta) * (1.0 - zeta) / denominator;
1691  functions[1][3] = (1.0 + eta) * (1.0 - zeta) / denominator;
1692  functions[1][4] = -(1.0 - eta) * (1.0 + zeta) / denominator;
1693  functions[1][5] = (1.0 - eta) * (1.0 + zeta) / denominator;
1694  functions[1][6] = -(1.0 + eta) * (1.0 + zeta) / denominator;
1695  functions[1][7] = (1.0 + eta) * (1.0 + zeta) / denominator;
1696 
1697  functions[2][0] = -(1.0 - xi) * (1.0 - zeta) / denominator;
1698  functions[2][1] = -(1.0 + xi) * (1.0 - zeta) / denominator;
1699  functions[2][2] = (1.0 - xi) * (1.0 - zeta) / denominator;
1700  functions[2][3] = (1.0 + xi) * (1.0 - zeta) / denominator;
1701  functions[2][4] = -(1.0 - xi) * (1.0 + zeta) / denominator;
1702  functions[2][5] = -(1.0 + xi) * (1.0 + zeta) / denominator;
1703  functions[2][6] = (1.0 - xi) * (1.0 + zeta) / denominator;
1704  functions[2][7] = (1.0 + xi) * (1.0 + zeta) / denominator;
1705 
1706  functions[3][0] = -(1.0 - xi) * (1.0 - eta) / denominator;
1707  functions[3][1] = -(1.0 + xi) * (1.0 - eta) / denominator;
1708  functions[3][2] = -(1.0 - xi) * (1.0 + eta) / denominator;
1709  functions[3][3] = -(1.0 + xi) * (1.0 + eta) / denominator;
1710  functions[3][4] = (1.0 - xi) * (1.0 - eta) / denominator;
1711  functions[3][5] = (1.0 + xi) * (1.0 - eta) / denominator;
1712  functions[3][6] = (1.0 - xi) * (1.0 + eta) / denominator;
1713  functions[3][7] = (1.0 + xi) * (1.0 + eta) / denominator;
1714 
1715 } // End GetInterpolationFunctions
1716 
1717 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1718 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort_permute(
1719  const typename Teuchos::Array<LocalOrdinal>::iterator& first1,
1720  const typename Teuchos::Array<LocalOrdinal>::iterator& last1,
1721  const typename Teuchos::Array<LocalOrdinal>::iterator& first2,
1722  const typename Teuchos::Array<LocalOrdinal>::iterator& /* last2 */) const {
1723  typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1724  DT n = last1 - first1;
1725  DT m = n / 2;
1726  DT z = Teuchos::OrdinalTraits<DT>::zero();
1727  while (m > z) {
1728  DT max = n - m;
1729  for (DT j = 0; j < max; j++) {
1730  for (DT k = j; k >= 0; k -= m) {
1731  if (first1[first2[k + m]] >= first1[first2[k]])
1732  break;
1733  std::swap(first2[k + m], first2[k]);
1734  }
1735  }
1736  m = m / 2;
1737  }
1738 } // End sh_sort_permute
1739 
1740 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1741 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort2(
1742  const typename Teuchos::Array<LocalOrdinal>::iterator& first1,
1743  const typename Teuchos::Array<LocalOrdinal>::iterator& last1,
1744  const typename Teuchos::Array<LocalOrdinal>::iterator& first2,
1745  const typename Teuchos::Array<LocalOrdinal>::iterator& /* last2 */) const {
1746  typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1747  DT n = last1 - first1;
1748  DT m = n / 2;
1749  DT z = Teuchos::OrdinalTraits<DT>::zero();
1750  while (m > z) {
1751  DT max = n - m;
1752  for (DT j = 0; j < max; j++) {
1753  for (DT k = j; k >= 0; k -= m) {
1754  if (first1[k + m] >= first1[k])
1755  break;
1756  std::swap(first1[k + m], first1[k]);
1757  std::swap(first2[k + m], first2[k]);
1758  }
1759  }
1760  m = m / 2;
1761  }
1762 } // End sh_sort2
1763 
1764 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1765 void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
1766  GetGIDLocalLexicographic(const GO i, const GO j, const GO k,
1767  const Array<LO> coarseNodeFineIndices, const RCP<GeometricData> myGeo,
1768  const LO myRankIndex, const LO pi, const LO pj, const LO /* pk */,
1769  const typename std::vector<std::vector<GO> >::iterator blockStart,
1770  const typename std::vector<std::vector<GO> >::iterator blockEnd,
1771  GO& myGID, LO& myPID, LO& myLID) const {
1772  LO ni = -1, nj = -1, li = -1, lj = -1, lk = -1;
1773  LO myRankGuess = myRankIndex;
1774  // We try to make a logical guess as to which PID owns the current coarse node
1775  if (i == 0 && myGeo->ghostInterface[0]) {
1776  --myRankGuess;
1777  } else if ((i == myGeo->ghostedCoarseNodesPerDir[0] - 1) && myGeo->ghostInterface[1]) {
1778  ++myRankGuess;
1779  }
1780  if (j == 0 && myGeo->ghostInterface[2]) {
1781  myRankGuess -= pi;
1782  } else if ((j == myGeo->ghostedCoarseNodesPerDir[1] - 1) && myGeo->ghostInterface[3]) {
1783  myRankGuess += pi;
1784  }
1785  if (k == 0 && myGeo->ghostInterface[4]) {
1786  myRankGuess -= pj * pi;
1787  } else if ((k == myGeo->ghostedCoarseNodesPerDir[2] - 1) && myGeo->ghostInterface[5]) {
1788  myRankGuess += pj * pi;
1789  }
1790  if (coarseNodeFineIndices[0] >= myGeo->meshData[myRankGuess][3] && coarseNodeFineIndices[0] <= myGeo->meshData[myRankGuess][4] && coarseNodeFineIndices[1] >= myGeo->meshData[myRankGuess][5] && coarseNodeFineIndices[1] <= myGeo->meshData[myRankGuess][6] && coarseNodeFineIndices[2] >= myGeo->meshData[myRankGuess][7] && coarseNodeFineIndices[2] <= myGeo->meshData[myRankGuess][8]) {
1791  myPID = myGeo->meshData[myRankGuess][0];
1792  ni = myGeo->meshData[myRankGuess][4] - myGeo->meshData[myRankGuess][3] + 1;
1793  nj = myGeo->meshData[myRankGuess][6] - myGeo->meshData[myRankGuess][5] + 1;
1794  li = coarseNodeFineIndices[0] - myGeo->meshData[myRankGuess][3];
1795  lj = coarseNodeFineIndices[1] - myGeo->meshData[myRankGuess][5];
1796  lk = coarseNodeFineIndices[2] - myGeo->meshData[myRankGuess][7];
1797  myLID = lk * nj * ni + lj * ni + li;
1798  myGID = myGeo->meshData[myRankGuess][9] + myLID;
1799  } else { // The guess failed, let us use the heavy artilery: std::find_if()
1800  // It could be interesting to monitor how many times this branch of the code gets
1801  // used as it is far more expensive than the above one...
1802  auto nodeRank = std::find_if(blockStart, blockEnd,
1803  [coarseNodeFineIndices](const std::vector<GO>& vec) {
1804  if (coarseNodeFineIndices[0] >= vec[3] && coarseNodeFineIndices[0] <= vec[4] && coarseNodeFineIndices[1] >= vec[5] && coarseNodeFineIndices[1] <= vec[6] && coarseNodeFineIndices[2] >= vec[7] && coarseNodeFineIndices[2] <= vec[8]) {
1805  return true;
1806  } else {
1807  return false;
1808  }
1809  });
1810  myPID = (*nodeRank)[0];
1811  ni = (*nodeRank)[4] - (*nodeRank)[3] + 1;
1812  nj = (*nodeRank)[6] - (*nodeRank)[5] + 1;
1813  li = coarseNodeFineIndices[0] - (*nodeRank)[3];
1814  lj = coarseNodeFineIndices[1] - (*nodeRank)[5];
1815  lk = coarseNodeFineIndices[2] - (*nodeRank)[7];
1816  myLID = lk * nj * ni + lj * ni + li;
1817  myGID = (*nodeRank)[9] + myLID;
1818  }
1819 } // End GetGIDLocalLexicographic
1820 
1821 } // namespace MueLu
1822 
1823 #define MUELU_GENERALGEOMETRICPFACTORY_SHORT
1824 #endif // MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
MueLu::GeneralGeometricPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const LO numDimension, const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1573
Teuchos_SerialDenseMatrix.hpp
MueLu::GeneralGeometricPFactory::GetCoarsePoints
void GetCoarsePoints(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:523
MueLu::GeneralGeometricPFactory::MakeGeneralGeometricP
void MakeGeneralGeometricP(RCP< GeometricData > myGeo, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &fCoords, const LO nnzP, const LO dofsPerNode, RCP< const Map > &stridedDomainMapP, RCP< Matrix > &Amat, RCP< Matrix > &P, RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &cCoords, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > coarseNodesGIDs, int interpolationOrder) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:719
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu::GeneralGeometricPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:107
Teuchos::Array::empty
bool empty() const
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:165
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LO
LocalOrdinal LO
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra_Matrix.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:15
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:34
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
Xpetra_IO.hpp
Teuchos_SerialDenseSolver.hpp
MueLu::GeneralGeometricPFactory::MeshLayoutInterface
void MeshLayoutInterface(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:288
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
Xpetra::MapFactory::Build
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Teuchos::SerialDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
MueLu::GeneralGeometricPFactory::sh_sort2
void sh_sort2(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1741
MueLu::GeneralGeometricPFactory::ComputeStencil
void ComputeStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, const int interpolationOrder, std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1530
MueLu::GeneralGeometricPFactory::GetGIDLocalLexicographic
void GetGIDLocalLexicographic(const GO i, const GO j, const GO k, const Array< LO > coarseNodeFineIndices, const RCP< GeometricData > myGeo, const LO myRankIndex, const LO pi, const LO pj, const LO pk, const typename std::vector< std::vector< GO > >::iterator blockStart, const typename std::vector< std::vector< GO > >::iterator blockEnd, GO &myGID, LO &myPID, LO &myLID) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1766
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::GeneralGeometricPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimension, const Teuchos::SerialDenseVector< LO, double > parameters, double functions[4][8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1661
MueLu_GeneralGeometricPFactory_decl.hpp
Xpetra::global_size_t
size_t global_size_t
Teuchos::Array::end
iterator end()
MueLu::GeneralGeometricPFactory::ComputeConstantInterpolationStencil
void ComputeConstantInterpolationStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1550
iter
TransListIter iter
MueLu::GeneralGeometricPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:69
MueLu::GeneralGeometricPFactory::sh_sort_permute
void sh_sort_permute(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1718
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::GeneralGeometricPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:101
MueLu_Monitor.hpp
Teuchos_SerialDenseVector.hpp
MueLu::GeneralGeometricPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:37
Xpetra_VectorFactory.hpp
NO
Node NO
Xpetra::StridedMapFactory::Build
static RCP< Xpetra::StridedMap< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:26
Teuchos::Array::iterator
std::vector< T >::iterator iterator
Teuchos::Array::begin
iterator begin()
MueLu::Level::DeclareInput
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
Definition: MueLu_Level.cpp:156
Xpetra_ImportFactory.hpp
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Definition: MueLu_Level.cpp:281
Teuchos::SerialDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)