MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_GeneralGeometricPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
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 
67 #include "MueLu_FactoryManagerBase.hpp"
68 #include "MueLu_GeneralGeometricPFactory_decl.hpp"
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>
76  RCP<const ParameterList> GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
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>
107  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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()),
118  Exceptions::RuntimeError,
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()),
131  Exceptions::RuntimeError,
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>
140  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel,
141  Level& coarseLevel) const {
142  return BuildP(fineLevel, coarseLevel);
143  }
144 
145  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
146  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& fineLevel,
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),
206  Exceptions::RuntimeError,
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),
211  Exceptions::RuntimeError,
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),
226  Exceptions::RuntimeError,
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),
244  Exceptions::RuntimeError,
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,
254  Exceptions::RuntimeError,
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>
331  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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>
566  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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>
777  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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, Xpetra::StaticProfile));
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>
1603  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeStencil(
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),
1610  Exceptions::RuntimeError,
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>
1623  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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>
1647  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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>
1736  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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>
1794  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort_permute(
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>
1821  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort2(
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>
1849  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, 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
MueLu::GeneralGeometricPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const LO numDimension, const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1648
Teuchos_SerialDenseMatrix.hpp
MueLu::GeneralGeometricPFactory::GetCoarsePoints
void GetCoarsePoints(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:567
MueLu::GeneralGeometricPFactory::MakeGeneralGeometricP
void MakeGeneralGeometricP(RCP< GeometricData > myGeo, const RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &fCoords, const LO nnzP, const LO dofsPerNode, RCP< const Map > &stridedDomainMapP, RCP< Matrix > &Amat, RCP< Matrix > &P, RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::coordinateType, LO, GO, NO > > &cCoords, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > coarseNodesGIDs, int interpolationOrder) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:778
Xpetra::StridedMapFactory::Build
static RCP< StridedMap > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, std::vector< size_t > &stridingInfo, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalOrdinal stridedBlockId=-1, GlobalOrdinal offset=0, LocalGlobal lg=Xpetra::GloballyDistributed)
MueLu_FactoryManagerBase.hpp
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu::GeneralGeometricPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:146
Teuchos::Array::empty
bool empty() const
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:209
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
LO
LocalOrdinal LO
Xpetra_Matrix.hpp
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu_MasterList.hpp
Teuchos::ArrayRCP::resize
void resize(const size_type n, const T &val=T())
Xpetra_IO.hpp
Teuchos_SerialDenseSolver.hpp
Xpetra_MultiVectorFactory.hpp
MueLu::GeneralGeometricPFactory::MeshLayoutInterface
void MeshLayoutInterface(const int interpolationOrder, const LO blkSize, RCP< const Map > fineCoordsMap, RCP< GeometricData > myGeometry, RCP< NodesIDs > ghostedCoarseNodes, Array< Array< GO > > &lCoarseNodesGIDs) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:332
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Teuchos::SerialDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
MueLu::GeneralGeometricPFactory::sh_sort2
void sh_sort2(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1821
MueLu::GeneralGeometricPFactory::ComputeStencil
void ComputeStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], const Array< Array< typename Teuchos::ScalarTraits< Scalar >::coordinateType > > coord, const int interpolationOrder, std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1603
MueLu::GeneralGeometricPFactory::GetGIDLocalLexicographic
void GetGIDLocalLexicographic(const GO i, const GO j, const GO k, const Array< LO > coarseNodeFineIndices, const RCP< GeometricData > myGeo, const LO myRankIndex, const LO pi, const LO pj, const LO pk, const typename std::vector< std::vector< GO > >::iterator blockStart, const typename std::vector< std::vector< GO > >::iterator blockEnd, GO &myGID, LO &myPID, LO &myLID) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1850
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::GeneralGeometricPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimension, const Teuchos::SerialDenseVector< LO, double > parameters, double functions[4][8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1737
MueLu_GeneralGeometricPFactory_decl.hpp
Xpetra::global_size_t
size_t global_size_t
Teuchos::Array::end
iterator end()
MueLu::GeneralGeometricPFactory::ComputeConstantInterpolationStencil
void ComputeConstantInterpolationStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], std::vector< double > &stencil) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1624
Xpetra_MapFactory.hpp
iter
TransListIter iter
MueLu::GeneralGeometricPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:108
MueLu::GeneralGeometricPFactory::sh_sort_permute
void sh_sort_permute(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1794
Teuchos::SerialDenseVector::size
int size(OrdinalType length_in)
Teuchos::SerialDenseSolver::setVectors
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
MueLu::GeneralGeometricPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:140
MueLu_Monitor.hpp
Teuchos_SerialDenseVector.hpp
MueLu::GeneralGeometricPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:76
Xpetra_VectorFactory.hpp
NO
Node NO
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
Xpetra_CrsMatrixWrap.hpp
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
Xpetra::MapFactory::Build
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int > > &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:63
Teuchos::Array::iterator
std::vector< T >::iterator iterator
Teuchos::Array::begin
iterator begin()
MueLu::Level::DeclareInput
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
Definition: MueLu_Level.cpp:150
Xpetra_ImportFactory.hpp
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Definition: MueLu_Level.hpp:332
Teuchos::SerialDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)