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