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