MueLu  Version of the Day
MueLu_GeneralGeometricPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
47 #define MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
48 
49 #include <stdlib.h>
50 #include <iomanip>
51 
52 
53 // #include <Teuchos_LAPACK.hpp>
54 #include <Teuchos_SerialDenseMatrix.hpp>
55 #include <Teuchos_SerialDenseVector.hpp>
56 #include <Teuchos_SerialDenseSolver.hpp>
57 
58 #include <Xpetra_CrsMatrixWrap.hpp>
59 #include <Xpetra_ImportFactory.hpp>
60 #include <Xpetra_Matrix.hpp>
61 #include <Xpetra_MapFactory.hpp>
62 #include <Xpetra_MultiVectorFactory.hpp>
63 #include <Xpetra_VectorFactory.hpp>
64 
65 #include <Xpetra_IO.hpp>
66 
67 #include "MueLu_FactoryManagerBase.hpp"
68 #include "MueLu_GeneralGeometricPFactory_decl.hpp"
69 
70 #include "MueLu_MasterList.hpp"
71 #include "MueLu_Monitor.hpp"
72 
73 namespace MueLu {
74 
75  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
76  RCP<const ParameterList> GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
77  RCP<ParameterList> validParamList = rcp(new ParameterList());
78 
79  // Coarsen can come in two forms, either a single char that will be interpreted as an integer which is used as the coarsening
80  // rate in every spatial dimentions, or it can be a longer string that will then be interpreted as an array of integers.
81  // By default coarsen is set as "{2}", hence a coarsening rate of 2 in every spatial dimension is the default setting!
82  validParamList->set<std::string > ("Coarsen", "{2}", "Coarsening rate in all spatial dimensions");
83  validParamList->set<int> ("order", 1, "Order of the interpolation scheme used");
84  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
85  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null, "Generating factory of the nullspace");
86  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null, "Generating factory for coorindates");
87  validParamList->set<RCP<const FactoryBase> >("gNodesPerDim", Teuchos::null, "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
88  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null, "Number of nodes per spatial dimmension provided by CoordinatesTransferFactory.");
89  validParamList->set<std::string >("axisPermutation", "{0,1,2}", "Assuming a global (x,y,z) orientation, local might be (z,y,x). This vector gives a permutation from global to local orientation.");
90 
91  return validParamList;
92  }
93 
94  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
95  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level& fineLevel, Level& coarseLevel) const {
96  Input(fineLevel, "A");
97  Input(fineLevel, "Nullspace");
98  Input(fineLevel, "Coordinates");
99  // Request the global number of nodes per dimensions
100  if(fineLevel.GetLevelID() == 0) {
101  if(fineLevel.IsAvailable("gNodesPerDim", NoFactory::get())) {
102  fineLevel.DeclareInput("gNodesPerDim", NoFactory::get(), this);
103  } else {
104  TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("gNodesPerDim", NoFactory::get()),
105  Exceptions::RuntimeError,
106  "gNodesPerDim was not provided by the user on level0!");
107  }
108  } else {
109  Input(fineLevel, "gNodesPerDim");
110  }
111 
112  // Request the local number of nodes per dimensions
113  if(fineLevel.GetLevelID() == 0) {
114  if(fineLevel.IsAvailable("lNodesPerDim", NoFactory::get())) {
115  fineLevel.DeclareInput("lNodesPerDim", NoFactory::get(), this);
116  } else {
117  TEUCHOS_TEST_FOR_EXCEPTION(fineLevel.IsAvailable("lNodesPerDim", NoFactory::get()),
118  Exceptions::RuntimeError,
119  "lNodesPerDim was not provided by the user on level0!");
120  }
121  } else {
122  Input(fineLevel, "lNodesPerDim");
123  }
124  }
125 
126  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
127  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level& fineLevel, Level& coarseLevel) const {
128  return BuildP(fineLevel, coarseLevel);
129  }
130 
131  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
132  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& fineLevel, Level& coarseLevel) const {
133  FactoryMonitor m(*this, "Build", coarseLevel);
134 
135  // Obtain general variables
136  using xdMV = Xpetra::MultiVector<double,LO,GO,NO>;
137  RCP<Matrix> A = Get< RCP<Matrix> > (fineLevel, "A");
138  RCP<MultiVector> fineNullspace = Get< RCP<MultiVector> > (fineLevel, "Nullspace");
139  RCP<xdMV> fineCoords = Get< RCP<xdMV> >(fineLevel, "Coordinates");
140  RCP<xdMV> coarseCoords;
141 
142  // Get user-provided coarsening rate parameter (constant over all levels)
143  const ParameterList& pL = GetParameterList();
144 
145  // collect general input data
146  LO blkSize = A->GetFixedBlockSize();
147  RCP<const Map> rowMap = A->getRowMap();
148  LO numDimensions = 0; // Number of spatial dimensions
149  Array<GO> gFineNodesPerDir(3); // Global number of fine points per direction
150  Array<GO> gCoarseNodesPerDir(3); // Global number of coarse points per direction
151  LO lNumFinePoints; // Local number of fine points
152  Array<LO> lFineNodesPerDir(3); // Local number of fine points per direction
153  Array<LO> lCoarseNodesPerDir(3); // Local number of coarse points per direction
154  Array<LO> mapDirL2G(3);
155  Array<LO> mapDirG2L(3);
156 
157  TEUCHOS_TEST_FOR_EXCEPTION(fineCoords==Teuchos::null, Exceptions::RuntimeError, "Coordinates cannot be accessed from fine level!");
158  numDimensions = fineCoords->getNumVectors();
159  lNumFinePoints = fineCoords->getLocalLength();
160 
161  // Get the number of points in each direction
162  if(fineLevel.GetLevelID() == 0) {
163  gFineNodesPerDir = fineLevel.Get<Array<GO> >("gNodesPerDim", NoFactory::get());
164  lFineNodesPerDir = fineLevel.Get<Array<LO> >("lNodesPerDim", NoFactory::get());
165  } else {
166  // Loading global number of nodes per diretions
167  gFineNodesPerDir = Get<Array<GO> >(fineLevel, "gNodesPerDim");
168 
169  // Loading local number of nodes per diretions
170  lFineNodesPerDir = Get<Array<LO> >(fineLevel, "lNodesPerDim");
171  }
172 
173  // Get the coarsening rate
174  std::string coarsenRate = pL.get<std::string>("Coarsen");
175  ArrayRCP<LO> coarseRate = arcp<LO>(3);
176  Teuchos::Array<LO> crates;
177  try {
178  crates = Teuchos::fromStringToArray<LO>(coarsenRate);
179  } catch(const Teuchos::InvalidArrayStringRepresentation e) {
180  GetOStream(Errors,-1) << " *** Coarsen must be a string convertible into an array! *** " << std::endl;
181  throw e;
182  }
183  TEUCHOS_TEST_FOR_EXCEPTION((crates.size()>1) && (crates.size()<numDimensions),
184  Exceptions::RuntimeError,
185  "Coarsen must have at least as many components as the number of spatial dimensions in the problem.");
186  for(LO i = 0; i < 3; ++i) {
187  if(i < numDimensions) {
188  if( crates.size()==1 ) {
189  coarseRate[i] = crates[0];
190  } else if( crates.size()==numDimensions ) {
191  coarseRate[i] = crates[i];
192  }
193  } else {
194  coarseRate[i] = 1;
195  }
196  }
197 
198  int interpolationOrder = pL.get<int>("order");
199  TEUCHOS_TEST_FOR_EXCEPTION((interpolationOrder < 0) || (interpolationOrder > 1),
200  Exceptions::RuntimeError,
201  "The interpolation order can only be set to 0 or 1.");
202 
203  // Get the axis permutation from Global axis to Local axis
204  std::string axisPermutation = pL.get<std::string>("axisPermutation");
205  try {
206  mapDirG2L = Teuchos::fromStringToArray<LO>(axisPermutation);
207  } catch(const Teuchos::InvalidArrayStringRepresentation e) {
208  GetOStream(Errors,-1) << " *** axisPermutation must be a string convertible into an array! *** " << std::endl;
209  throw e;
210  }
211  for(LO i = 0; i < numDimensions; ++i) {
212  TEUCHOS_TEST_FOR_EXCEPTION(mapDirG2L[i] > numDimensions,
213  Exceptions::RuntimeError,
214  "axis permutation values must all be less than the number of spatial dimensions.");
215  mapDirL2G[mapDirG2L[i]] = i;
216  }
217 
218  // Find the offsets for the nodes on the current processor, this tells us the position of the first node on the partition compare to the
219  // closest lower coarse point. This information is needed to know who are the coarse nodes that need to be used for interpolation and by
220  // extension this tells us what are the ghostnodes we need for the computation.
221  //
222  // x - o - o - x
223  // | | | |
224  // o - o - o - o
225  // | | | |
226  // o - o - * - o
227  // | | | |
228  // x - o - o - x
229  //
230  // x are coarse points, o are fine points and * is the first point on the local processor.
231  // the offsets of * are 2 and 1
232  GO minGlobalIndex, maxGlobalIndex, rem;
233  GO gIndices[6] = {0, 0, 0, 0, 0, 0};
234  LO offsets[6] = {0, 0, 0, 0, 0, 0};
235  RCP<const Map> fineCoordsMap = fineCoords->getMap();
236  minGlobalIndex = fineCoordsMap->getMinGlobalIndex();
237  maxGlobalIndex = fineCoordsMap->getMaxGlobalIndex();
238  if(numDimensions == 1) {
239  gIndices[0] = minGlobalIndex;
240  offsets[0] = Teuchos::as<LO>(gIndices[0]) % coarseRate[0];
241 
242  gIndices[3] = maxGlobalIndex;
243  offsets[3] = Teuchos::as<LO>(gIndices[3]) % coarseRate[0];
244  } else if(numDimensions == 2) {
245  gIndices[1] = minGlobalIndex / Teuchos::as<GO>(gFineNodesPerDir[0]);
246  offsets[1] = Teuchos::as<LO>(gIndices[1]) % coarseRate[1];
247  gIndices[0] = minGlobalIndex % Teuchos::as<GO>(gFineNodesPerDir[0]);
248  offsets[0] = Teuchos::as<LO>(gIndices[0]) % coarseRate[0];
249 
250  gIndices[4] = maxGlobalIndex / Teuchos::as<GO>(gFineNodesPerDir[0]);
251  offsets[4] = Teuchos::as<LO>(gIndices[4]) % coarseRate[1];
252  gIndices[3] = maxGlobalIndex % Teuchos::as<GO>(gFineNodesPerDir[0]);
253  offsets[3] = Teuchos::as<LO>(gIndices[3]) % coarseRate[3];
254  } else if(numDimensions == 3) {
255  gIndices[2] = minGlobalIndex / Teuchos::as<GO>(gFineNodesPerDir[0]*gFineNodesPerDir[1]);
256  rem = minGlobalIndex % Teuchos::as<GO>(gFineNodesPerDir[0]*gFineNodesPerDir[1]);
257  offsets[2] = Teuchos::as<LO>(gIndices[2]) % coarseRate[2];
258  gIndices[1] = rem / Teuchos::as<GO>(gFineNodesPerDir[0]);
259  offsets[1] = Teuchos::as<LO>(gIndices[1]) % coarseRate[1];
260  gIndices[0] = rem % Teuchos::as<GO>(gFineNodesPerDir[0]);
261  offsets[0] = Teuchos::as<LO>(gIndices[0]) % coarseRate[0];
262 
263  gIndices[5] = maxGlobalIndex / Teuchos::as<GO>(gFineNodesPerDir[0]*gFineNodesPerDir[1]);
264  rem = maxGlobalIndex % Teuchos::as<GO>(gFineNodesPerDir[0]*gFineNodesPerDir[1]);
265  offsets[5] = Teuchos::as<LO>(gIndices[5]) % coarseRate[2];
266  gIndices[4] = rem / Teuchos::as<GO>(gFineNodesPerDir[0]);
267  offsets[4] = Teuchos::as<LO>(gIndices[4]) % coarseRate[1];
268  gIndices[3] = rem % Teuchos::as<GO>(gFineNodesPerDir[0]);
269  offsets[3] = Teuchos::as<LO>(gIndices[3]) % coarseRate[0];
270  }
271  /* At this point the local geometrical discovery is over */
272 
273  // Check if the partition contains nodes on a boundary, if so that boundary (face, line or point) will not require ghost nodes.
274  // Always include a layer if ghost nodes for inner interface since even for faces with coarse nodes, an outward curvature might
275  // require ghost nodes to compute a proper interpolation operator. Curvature could be check localy to avoid including extra
276  // ghost nodes...
277  bool ghostInterface[6] = {false, false, false, false, false, false};
278  for(LO i=0; i < numDimensions; ++i) {
279  if( gIndices[i] != 0) {
280  ghostInterface[2*i]=true;
281  }
282  if( (gIndices[i]+lFineNodesPerDir[mapDirL2G[i]]) != gFineNodesPerDir[i] ) {
283  ghostInterface[2*i+1]=true;
284  }
285  }
286 
287  /* Here one element can represent either the degenerate case of one node or the more general case of two nodes.
288  i.e. x---x is a 1D element with two nodes and x is a 1D element with one node. This helps generating a 3D space from tensorial products...
289  A good way to handle this would be to generalize the algorithm to take into account the discretization order used in each direction,
290  at least in the FEM sense, since a 0 degree discretization will have a unique node per element. This way a 1D discretization can be viewed
291  as a 3D problem with one 0 degree element in the y direction and one 0 degre element in the z direction. */
292  LO endRate[3] = {0, 0, 0};
293  /* /!\ Operations below are aftecting both local and global values that have two different orientations. /!\
294  orientations can be interchanged using mapDirG2L and mapDirL2G. coarseRate, endRate and offsets
295  are in the global basis, as well as all the variables starting with a g, while the variables
296  /!\ starting with an l are in the local basis. /!\ */
297  for(LO i = 0; i < 3; ++i) {
298  if(i < numDimensions) {
299  // This array is passed to the RAPFactory and eventually becomes gFineNodePerDir on the next level.
300  gCoarseNodesPerDir[i] = (gFineNodesPerDir[i] - 1) / coarseRate[i];
301  endRate[i] = Teuchos::as<LO>((gFineNodesPerDir[i] - 1) % coarseRate[i]);
302  if(endRate[i] == 0) {
303  endRate[i] = coarseRate[i];
304  ++gCoarseNodesPerDir[i];
305  } else {
306  gCoarseNodesPerDir[i] += 2;
307  }
308  } else {
309  endRate[i] = 1;
310  gCoarseNodesPerDir[i] = 1;
311  }
312  }
313 
314  for(LO i = 0; i < 3; ++i) {
315  if(i < numDimensions) {
316  // Check whether the partition includes the "end" of the mesh which means that endRate will apply.
317  // Also make sure that endRate is not 0 which means that the mesh does not require a particular treatment
318  // at the boundaries.
319  if( (gIndices[mapDirG2L[i]] + lFineNodesPerDir[i]) == gFineNodesPerDir[mapDirG2L[i]] ) {
320  lCoarseNodesPerDir[i] = (lFineNodesPerDir[i] - endRate[mapDirG2L[i]] + offsets[mapDirG2L[i]] - 1) / coarseRate[mapDirG2L[i]] + 1;
321  if(offsets[mapDirG2L[i]] == 0) {++lCoarseNodesPerDir[i];}
322  } else {
323  lCoarseNodesPerDir[i] = (lFineNodesPerDir[i] + offsets[mapDirG2L[i]] - 1) / coarseRate[mapDirG2L[i]];
324  if(offsets[mapDirG2L[i]] == 0) {++lCoarseNodesPerDir[i];}
325  }
326  } else {
327  lCoarseNodesPerDir[i] = 1;
328  }
329  if(lFineNodesPerDir[i] < 1) {lCoarseNodesPerDir[i] = 0;}
330  }
331 
332  // Assuming linear interpolation, each fine point has contribution from 8 coarse points
333  // and each coarse point value gets injected.
334  // For systems of PDEs we assume that all dofs have the same P operator.
335  LO lNumCoarsePoints = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];
336  LO nTerms = 8*lNumFinePoints - 7*lNumCoarsePoints;
337  nTerms=nTerms*blkSize;
338 
339  // For each direction, determine how many ghost points are required.
340  LO numGhosts = 0;
341  const LO complementaryIndices[3][2] = {{1,2}, {0,2}, {0,1}};
342  for(LO i = 0; i < 3; ++i) {
343  LO tmp = 0;
344  // Check whether a face in direction i needs ghost nodes
345  if(ghostInterface[2*i] || ghostInterface[2*i+1]) {
346  if(i == 0) {tmp = lCoarseNodesPerDir[mapDirL2G[1]]*lCoarseNodesPerDir[mapDirL2G[2]];}
347  if(i == 1) {tmp = lCoarseNodesPerDir[mapDirL2G[0]]*lCoarseNodesPerDir[mapDirL2G[2]];}
348  if(i == 2) {tmp = lCoarseNodesPerDir[mapDirL2G[0]]*lCoarseNodesPerDir[mapDirL2G[1]];}
349  }
350  // If both faces in direction i need nodes, double the number of ghost nodes
351  if(ghostInterface[2*i] && ghostInterface[2*i+1]) {tmp = 2*tmp;}
352  numGhosts += tmp;
353 
354  // The corners and edges need to be checked in 2D / 3D to add more ghosts nodes
355  for(LO j = 0; j < 2; ++j) {
356  for(LO k = 0; k < 2; ++k) {
357  // Check if two adjoining faces need ghost nodes and then add their common edge
358  if(ghostInterface[2*complementaryIndices[i][0]+j] && ghostInterface[2*complementaryIndices[i][1]+k]) {
359  numGhosts += lCoarseNodesPerDir[mapDirL2G[i]];
360  // Add corners if three adjoining faces need ghost nodes,
361  // but add them only once! Hence when i == 0.
362  if(ghostInterface[2*i] && (i == 0)) { numGhosts += 1; }
363  if(ghostInterface[2*i+1] && (i == 0)) { numGhosts += 1; }
364  }
365  }
366  }
367  }
368 
369  // Build a list of GIDs to import the required ghost nodes.
370  // The ordering of the ghosts nodes will be as natural as possible,
371  // i.e. it should follow the ordering of the mesh.
372  //
373  // Saddly we have to more or less redo what was just done to figure out the value of numGhosts,
374  // there might be some optimization possibility here...
375  Array<GO> ghostsGIDs(numGhosts);
376  LO countGhosts = 0;
377  // Get the GID of the first point on the current partition.
378  GO startingGID = minGlobalIndex;
379  Array<GO> startingIndices(3);
380  // We still want ghost nodes even if have with a 0 offset,
381  // except when we are on a boundary
382  if(ghostInterface[4] && (offsets[2] == 0)) {
383  if(gIndices[2] + coarseRate[2] > gFineNodesPerDir[2]) {
384  startingGID -= endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
385  } else {
386  startingGID -= coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
387  }
388  } else {
389  startingGID -= offsets[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
390  }
391  if(ghostInterface[2] && (offsets[1] == 0)) {
392  if(gIndices[1] + coarseRate[1] > gFineNodesPerDir[1]) {
393  startingGID -= endRate[1]*gFineNodesPerDir[0];
394  } else {
395  startingGID -= coarseRate[1]*gFineNodesPerDir[0];
396  }
397  } else {
398  startingGID -= offsets[1]*gFineNodesPerDir[0];
399  }
400  if(ghostInterface[0] && (offsets[0] == 0)) {
401  if(gIndices[0] + coarseRate[0] > gFineNodesPerDir[0]) {
402  startingGID -= endRate[0];
403  } else {
404  startingGID -= coarseRate[0];
405  }
406  } else {
407  startingGID -= offsets[0];
408  }
409 
410  { // scope for tmp
411  GO tmp;
412  startingIndices[2] = startingGID / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
413  tmp = startingGID % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
414  startingIndices[1] = tmp / gFineNodesPerDir[0];
415  startingIndices[0] = tmp % gFineNodesPerDir[0];
416  }
417 
418  GO ghostOffset[3] = {0, 0, 0};
419  LO lengthZero = lCoarseNodesPerDir[mapDirL2G[0]];
420  LO lengthOne = lCoarseNodesPerDir[mapDirL2G[1]];
421  LO lengthTwo = lCoarseNodesPerDir[mapDirL2G[2]];
422  if(ghostInterface[0]) {++lengthZero;}
423  if(ghostInterface[1]) {++lengthZero;}
424  if(ghostInterface[2]) {++lengthOne;}
425  if(ghostInterface[3]) {++lengthOne;}
426  if(ghostInterface[4]) {++lengthTwo;}
427  if(ghostInterface[5]) {++lengthTwo;}
428 
429  // First check the bottom face as it will have the lowest GIDs
430  if(ghostInterface[4]) {
431  ghostOffset[2] = startingGID;
432  for(LO j = 0; j < lengthOne; ++j) {
433  if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
434  ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
435  } else {
436  ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
437  }
438  for(LO k = 0; k < lengthZero; ++k) {
439  if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
440  ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
441  } else {
442  ghostOffset[0] = k*coarseRate[0];
443  }
444  // If the partition includes a changed rate at the edge, ghost nodes need to be picked carefully.
445  // This if statement is repeated each time a ghost node is selected.
446  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
447  ++countGhosts;
448  }
449  }
450  }
451 
452  // Sweep over the lCoarseNodesPerDir[2] coarse layers in direction 2 and gather necessary ghost nodes
453  // located on these layers.
454  for(LO i = 0; i < lengthTwo; ++i) {
455  // Exclude the cases where ghost nodes exists on the faces in directions 2, these faces are swept
456  // seperatly for efficiency.
457  if( !((i == lengthTwo-1) && ghostInterface[5]) && !((i == 0) && ghostInterface[4]) ) {
458  // Set the ghostOffset in direction 2 taking into account a possible endRate different
459  // from the regular coarseRate.
460  if( (i == lengthTwo-1) && !ghostInterface[5] ) {
461  ghostOffset[2] = startingGID + ((i-1)*coarseRate[2] + endRate[2])*gFineNodesPerDir[1]*gFineNodesPerDir[0];
462  } else {
463  ghostOffset[2] = startingGID + i*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
464  }
465  for(LO j = 0; j < lengthOne; ++j) {
466  if( (j == 0) && ghostInterface[2] ) {
467  for(LO k = 0; k < lengthZero; ++k) {
468  if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
469  if(k == 0) {
470  ghostOffset[0] = endRate[0];
471  } else {
472  ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
473  }
474  } else {
475  ghostOffset[0] = k*coarseRate[0];
476  }
477  // In this case j == 0 so ghostOffset[1] is zero and can be ignored
478  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[0];
479  ++countGhosts;
480  }
481  } else if( (j == lengthOne-1) && ghostInterface[3] ) {
482  // Set the ghostOffset in direction 1 taking into account a possible endRate different
483  // from the regular coarseRate.
484  if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
485  ghostOffset[1] = ((j-1)*coarseRate[1] + endRate[1])*gFineNodesPerDir[0];
486  } else {
487  ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
488  }
489  for(LO k = 0; k < lengthZero; ++k) {
490  if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {
491  ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
492  } else {
493  ghostOffset[0] = k*coarseRate[0];
494  }
495  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
496  ++countGhosts;
497  }
498  } else {
499  // Set ghostOffset[1] for side faces sweep
500  if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) {
501  ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
502  } else {
503  ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
504  }
505 
506  // Set ghostOffset[0], ghostsGIDs and countGhosts
507  if(ghostInterface[0]) { // In that case ghostOffset[0]==0, so we can ignore it
508  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1];
509  ++countGhosts;
510  }
511  if(ghostInterface[1]) { // Grab ghost point at the end of direction 0.
512  if( (startingIndices[0] + (lengthZero-1)*coarseRate[0]) > gFineNodesPerDir[0] - 1 ) {
513  if(lengthZero > 2) {
514  ghostOffset[0] = (lengthZero-2)*coarseRate[0] + endRate[0];
515  } else {
516  ghostOffset[0] = endRate[0];
517  }
518  } else {
519  ghostOffset[0] = (lengthZero-1)*coarseRate[0];
520  }
521  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
522  ++countGhosts;
523  }
524  }
525  }
526  }
527  }
528 
529  // Finally check the top face
530  if(ghostInterface[5]) {
531  if( startingIndices[2] + (lengthTwo-1)*coarseRate[2] + 1 > gFineNodesPerDir[2] ) {
532  ghostOffset[2] = startingGID + ((lengthTwo-2)*coarseRate[2] + endRate[2])*gFineNodesPerDir[1]*gFineNodesPerDir[0];
533  } else {
534  ghostOffset[2] = startingGID + (lengthTwo-1)*coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
535  }
536  for(LO j = 0; j < lengthOne; ++j) {
537  if( (j == lengthOne-1) && (startingIndices[1] + j*coarseRate[1] + 1 > gFineNodesPerDir[1]) ) { // && !ghostInterface[3] ) {
538  ghostOffset[1] = ( (j-1)*coarseRate[1] + endRate[1] )*gFineNodesPerDir[0];
539  } else {
540  ghostOffset[1] = j*coarseRate[1]*gFineNodesPerDir[0];
541  }
542  for(LO k = 0; k < lengthZero; ++k) {
543  if( (k == lengthZero-1) && (startingIndices[0] + k*coarseRate[0] + 1 > gFineNodesPerDir[0]) ) {// && !ghostInterface[1] ) {
544  ghostOffset[0] = (k-1)*coarseRate[0] + endRate[0];
545  } else {
546  ghostOffset[0] = k*coarseRate[0];
547  }
548  ghostsGIDs[countGhosts] = ghostOffset[2] + ghostOffset[1] + ghostOffset[0];
549  ++countGhosts;
550  }
551  }
552  }
553 
554  // All that is left to do is loop over NCpts and:
555  // - extract coarse points coordiate for coarseCoords
556  // - get coordinates for current stencil computation
557  // - compute current stencil
558  // - compute row and column indices for stencil entries
559  RCP<const Map> stridedDomainMapP;
560  RCP<Matrix> P;
561  MakeGeneralGeometricP(numDimensions, mapDirL2G, mapDirG2L, lFineNodesPerDir, lCoarseNodesPerDir, gCoarseNodesPerDir,
562  gFineNodesPerDir, coarseRate, endRate, offsets, ghostInterface, fineCoords, nTerms, blkSize,
563  stridedDomainMapP, A, P, coarseCoords, ghostsGIDs, interpolationOrder);
564 
565  // set StridingInformation of P
566  if (A->IsView("stridedMaps") == true) {
567  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMapP);
568  } else {
569  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMapP);
570  }
571 
572  // store the transfer operator and node coordinates on coarse level
573  Set(coarseLevel, "P", P);
574  Set(coarseLevel, "coarseCoordinates", coarseCoords);
575  Set<Array<GO> >(coarseLevel, "gCoarseNodesPerDim", gCoarseNodesPerDir);
576  Set<Array<LO> >(coarseLevel, "lCoarseNodesPerDim", lCoarseNodesPerDir);
577 
578  // rst: null space might get scaled here ... do we care. We could just inject at the cpoints, but I don't
579  // feel that this is needed.
580  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(), fineNullspace->getNumVectors());
581  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(), Teuchos::ScalarTraits<SC>::zero());
582  Set(coarseLevel, "Nullspace", coarseNullspace);
583 
584  }
585 
586  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
587  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::MakeGeneralGeometricP(LO const ndm, const Array<LO> mapDirL2G, const Array<LO> mapDirG2L, const Array<LO> lFineNodesPerDir, const Array<LO> lCoarseNodesPerDir, Array<GO> gCoarseNodesPerDir, Array<GO> gFineNodesPerDir, ArrayRCP<LO> const coarseRate, LO const endRate[3], LO const offsets[6], bool const ghostInterface[6], const RCP<Xpetra::MultiVector<double,LO,GO,Node> >& fineCoords, LO const nnzP, LO const dofsPerNode, RCP<const Map>& stridedDomainMapP, RCP<Matrix> & Amat, RCP<Matrix>& P, RCP<Xpetra::MultiVector<double,LO,GO,Node> >& coarseCoords, Array<GO> ghostsGIDs, int interpolationOrder) const {
588 
589  /*
590  *
591  *
592  *
593  *
594  *
595  *
596  *
597  *
598  *
599  *
600  *
601  * On termination, return the number of local prolongator columns owned by
602  * this processor.
603  *
604  * Input
605  * =====
606  * nNodes Number of fine level Blk Rows owned by this processor
607  * coarseRate Rate of coarsening in each spatial direction.
608  * endRate Rate of coarsening in each spatial direction for the last
609  * nodes in the mesh where an adaptive coarsening rate is
610  * required.
611  * nTerms Number of nonzero entries in the prolongation matrix.
612  * dofsPerNode Number of degrees-of-freedom per mesh node.
613  *
614  * Output
615  * =====
616  * So far nothing...
617  */
618 
619  using xdMV = Xpetra::MultiVector<double,LO,GO,NO>;
620 
621  LO lNumFineNodes = lFineNodesPerDir[0]*lFineNodesPerDir[1]*lFineNodesPerDir[2];
622  LO lNumCoarseNodes = lCoarseNodesPerDir[0]*lCoarseNodesPerDir[1]*lCoarseNodesPerDir[2];
623  GO gNumCoarseNodes = gCoarseNodesPerDir[0]*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[2];
624  LO lNumGhostNodes = ghostsGIDs.size();
625  GO numGloCols = dofsPerNode*gNumCoarseNodes;
626 
627  // Build the required column map for the prolongator operator.
628  // This requies to find the GIDs of the coarse nodes on the coarse mesh,
629  // including the ghost nodes, on the local partition.
630  // Note: the ordering of the coarse nodes is assumed
631  // to be the same as the ordering of the fine nodes.
632 
633  GO gStartIndices[3];
634  RCP<const Map> fineCoordsMap = fineCoords->getMap();
635  // Compute the global indices of the first node on the partition.
636  { // Scope for dummy
637  gStartIndices[2] = fineCoordsMap->getMinGlobalIndex() / (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
638  GO dummy = fineCoordsMap->getMinGlobalIndex() % (gFineNodesPerDir[1]*gFineNodesPerDir[0]);
639  gStartIndices[1] = dummy / gFineNodesPerDir[0];
640  gStartIndices[0] = dummy % gFineNodesPerDir[0];
641  }
642 
643  Array<GO> gCoarseNodesGIDs(lNumCoarseNodes);
644  LO currentNode, offset2, offset1, offset0;
645  // Find the GIDs of the coarse nodes on the partition.
646  for(LO ind2 = 0; ind2 < lCoarseNodesPerDir[mapDirL2G[2]]; ++ind2) {
647  if(offsets[2] == 0) {
648  offset2 = gStartIndices[2];
649  } else {
650  if(gStartIndices[2] + endRate[2] - offsets[2] == gFineNodesPerDir[2] - 1) {
651  offset2 = gStartIndices[2] + endRate[2] - offsets[2];
652  } else {
653  offset2 = gStartIndices[2] + coarseRate[2] - offsets[2];
654  }
655  }
656  if(offset2 + ind2*coarseRate[2] > gFineNodesPerDir[2] - 1) {
657  offset2 += (ind2 - 1)*coarseRate[2] + endRate[2];
658  } else {
659  offset2 += ind2*coarseRate[2];
660  }
661  offset2 = offset2*gFineNodesPerDir[1]*gFineNodesPerDir[0];
662 
663  for(LO ind1 = 0; ind1 < lCoarseNodesPerDir[mapDirL2G[1]]; ++ind1) {
664  if(offsets[1] == 0) {
665  offset1 = gStartIndices[1];
666  } else {
667  if(gStartIndices[1] + endRate[1] - offsets[1] == gFineNodesPerDir[1] - 1) {
668  offset1 = gStartIndices[1] + endRate[1] - offsets[1];
669  } else {
670  offset1 = gStartIndices[1] + coarseRate[1] - offsets[1];
671  }
672  }
673  if(offset1 + ind1*coarseRate[1] > gFineNodesPerDir[1] - 1) {
674  offset1 += (ind1 - 1)*coarseRate[1] + endRate[1];
675  } else {
676  offset1 += ind1*coarseRate[1];
677  }
678  offset1 = offset1*gFineNodesPerDir[0];
679  for(LO ind0 = 0; ind0 < lCoarseNodesPerDir[mapDirL2G[0]]; ++ind0) {
680  offset0 = gStartIndices[0] - offsets[0];
681  if(offsets[0] == 0) {
682  offset0 += ind0*coarseRate[0];
683  } else {
684  offset0 += (ind0 + 1)*coarseRate[0];
685  }
686  if(offset0 > gFineNodesPerDir[0] - 1) {offset0 += endRate[0] - coarseRate[0];}
687 
688  currentNode = ind2*lCoarseNodesPerDir[mapDirL2G[1]]*lCoarseNodesPerDir[mapDirL2G[0]]
689  + ind1*lCoarseNodesPerDir[mapDirL2G[0]]
690  + ind0;
691  gCoarseNodesGIDs[currentNode] = offset2 + offset1 + offset0;
692  }
693  }
694  }
695 
696  // Actual loop over all the coarse/ghost nodes to find their index on the coarse mesh
697  // and the corresponding dofs that will need to be added to colMapP.
698  Array<GO> colGIDs(dofsPerNode*(lNumCoarseNodes+lNumGhostNodes));
699  Array<GO> coarseNodesGIDs(lNumCoarseNodes);
700  GO fineNodesPerCoarseSlab = coarseRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
701  GO fineNodesEndCoarseSlab = endRate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0];
702  GO fineNodesPerCoarsePlane = coarseRate[1]*gFineNodesPerDir[0];
703  GO fineNodesEndCoarsePlane = endRate[1]*gFineNodesPerDir[0];
704  GO coarseNodesPerCoarseLayer = gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0];
705  GO gCoarseNodeOnCoarseGridGID;
706  LO gInd[3], lCol;
707  Array<int> ghostsPIDs (lNumGhostNodes);
708  Array<LO> ghostsLIDs (lNumGhostNodes);
709  Array<LO> ghostsPermut(lNumGhostNodes);
710  for(LO k = 0; k < lNumGhostNodes; ++k) {ghostsPermut[k] = k;}
711  fineCoordsMap->getRemoteIndexList(ghostsGIDs, ghostsPIDs, ghostsLIDs);
712  sh_sort_permute(ghostsPIDs.begin(),ghostsPIDs.end(), ghostsPermut.begin(),ghostsPermut.end());
713 
714  { // scope for tmpInds, tmpVars and tmp.
715  GO tmpInds[3], tmpVars[2];
716  LO tmp;
717  // Loop over the coarse nodes of the partition and add them to colGIDs
718  // that will be used to construct the column and domain maps of P as well
719  // as to construct the coarse coordinates map.
720  for(LO col = 0; col < lNumCoarseNodes; ++col) {
721  if((endRate[2] != coarseRate[2]) && (gCoarseNodesGIDs[col] > (gCoarseNodesPerDir[2] - 2)*fineNodesPerCoarseSlab + fineNodesEndCoarseSlab - 1)) {
722  tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab + 1;
723  tmpVars[0] = gCoarseNodesGIDs[col] - (tmpInds[2] - 1)*fineNodesPerCoarseSlab - fineNodesEndCoarseSlab;
724  } else {
725  tmpInds[2] = gCoarseNodesGIDs[col] / fineNodesPerCoarseSlab;
726  tmpVars[0] = gCoarseNodesGIDs[col] % fineNodesPerCoarseSlab;
727  }
728  if((endRate[1] != coarseRate[1]) && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane + fineNodesEndCoarsePlane - 1)) {
729  tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
730  tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane - fineNodesEndCoarsePlane;
731  } else {
732  tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
733  tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
734  }
735  if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
736  tmpInds[0] = gCoarseNodesPerDir[0] - 1;
737  } else {
738  tmpInds[0] = tmpVars[1] / coarseRate[0];
739  }
740  gInd[2] = col / (lCoarseNodesPerDir[mapDirG2L[1]]*lCoarseNodesPerDir[mapDirG2L[0]]);
741  tmp = col % (lCoarseNodesPerDir[mapDirG2L[1]]*lCoarseNodesPerDir[mapDirG2L[0]]);
742  gInd[1] = tmp / lCoarseNodesPerDir[mapDirG2L[0]];
743  gInd[0] = tmp % lCoarseNodesPerDir[mapDirG2L[0]];
744  lCol = gInd[mapDirG2L[2]]*(lCoarseNodesPerDir[1]*lCoarseNodesPerDir[0]) + gInd[mapDirG2L[1]]*lCoarseNodesPerDir[0] + gInd[mapDirG2L[0]];
745  gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
746  coarseNodesGIDs[lCol] = gCoarseNodeOnCoarseGridGID;
747  for(LO dof = 0; dof < dofsPerNode; ++dof) {
748  colGIDs[dofsPerNode*lCol + dof] = dofsPerNode*gCoarseNodeOnCoarseGridGID + dof;
749  }
750  }
751  // Now loop over the ghost nodes of the partition to add them to colGIDs
752  // since they will need to be included in the column map of P
753  for(LO col = lNumCoarseNodes; col < lNumCoarseNodes + lNumGhostNodes; ++col) {
754  if((endRate[2] != coarseRate[2]) && (ghostsGIDs[ghostsPermut[col - lNumCoarseNodes]] > (gCoarseNodesPerDir[2] - 2)*fineNodesPerCoarseSlab + fineNodesEndCoarseSlab - 1)) {
755  tmpInds[2] = ghostsGIDs[ghostsPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab + 1;
756  tmpVars[0] = ghostsGIDs[ghostsPermut[col - lNumCoarseNodes]] - (tmpInds[2] - 1)*fineNodesPerCoarseSlab - fineNodesEndCoarseSlab;
757  } else {
758  tmpInds[2] = ghostsGIDs[ghostsPermut[col - lNumCoarseNodes]] / fineNodesPerCoarseSlab;
759  tmpVars[0] = ghostsGIDs[ghostsPermut[col - lNumCoarseNodes]] % fineNodesPerCoarseSlab;
760  }
761  if((endRate[1] != coarseRate[1]) && (tmpVars[0] > (gCoarseNodesPerDir[1] - 2)*fineNodesPerCoarsePlane + fineNodesEndCoarsePlane - 1)) {
762  tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane + 1;
763  tmpVars[1] = tmpVars[0] - (tmpInds[1] - 1)*fineNodesPerCoarsePlane - fineNodesEndCoarsePlane;
764  } else {
765  tmpInds[1] = tmpVars[0] / fineNodesPerCoarsePlane;
766  tmpVars[1] = tmpVars[0] % fineNodesPerCoarsePlane;
767  }
768  if(tmpVars[1] == gFineNodesPerDir[0] - 1) {
769  tmpInds[0] = gCoarseNodesPerDir[0] - 1;
770  } else {
771  tmpInds[0] = tmpVars[1] / coarseRate[0];
772  }
773  gCoarseNodeOnCoarseGridGID = tmpInds[2]*coarseNodesPerCoarseLayer + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
774  for(LO dof = 0; dof < dofsPerNode; ++dof) {
775  colGIDs[dofsPerNode*col + dof] = dofsPerNode*gCoarseNodeOnCoarseGridGID + dof;
776  }
777  }
778  } // End of scope for tmpInds, tmpVars and tmp
779 
780  // Build maps necessary to create and fill complete the prolongator
781  // note: rowMapP == rangeMapP and colMapP != domainMapP
782  RCP<const Map> rowMapP = Amat->getDomainMap();
783 
784  RCP<const Map> domainMapP = Xpetra::MapFactory<LO,GO,NO>::Build(rowMapP->lib(),
785  numGloCols,
786  colGIDs.view(0, dofsPerNode*lNumCoarseNodes),
787  rowMapP->getIndexBase(),
788  rowMapP->getComm());
789 
790  RCP<const Map> colMapP = Xpetra::MapFactory<LO,GO,NO>::Build(rowMapP->lib(),
791  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
792  colGIDs.view(0, colGIDs.size()),
793  rowMapP->getIndexBase(),
794  rowMapP->getComm());
795 
796  std::vector<size_t> strideInfo(1);
797  strideInfo[0] = dofsPerNode;
798  stridedDomainMapP = Xpetra::StridedMapFactory<LO,GO,NO>::Build(domainMapP, strideInfo);
799 
800  // Build the map for the coarse level coordinates
801  RCP<const Map> coarseCoordsMap = MapFactory::Build (fineCoordsMap->lib(),
802  gNumCoarseNodes,
803  coarseNodesGIDs.view(0, lNumCoarseNodes),
804  fineCoordsMap->getIndexBase(),
805  rowMapP->getComm());
806 
807  // Do the actual import using the fineCoordsMap
808  RCP<const Map> ghostMap = Xpetra::MapFactory<LO,GO,NO>::Build(fineCoordsMap->lib(), Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(), ghostsGIDs.view(0, lNumGhostNodes),
809  fineCoordsMap->getIndexBase(), rowMapP->getComm());
810  RCP<const Import> importer = ImportFactory::Build(fineCoordsMap, ghostMap);
811  RCP<xdMV> ghostCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(ghostMap, ndm);
812  ghostCoords->doImport(*fineCoords, *importer, Xpetra::INSERT);
813 
814  P = rcp(new CrsMatrixWrap(rowMapP, colMapP, 0, Xpetra::StaticProfile));
815  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
816 
817  ArrayRCP<size_t> iaP;
818  ArrayRCP<LO> jaP;
819  ArrayRCP<SC> valP;
820 
821  PCrs->allocateAllValues(nnzP, iaP, jaP, valP);
822 
823  ArrayView<size_t> ia = iaP();
824  ArrayView<LO> ja = jaP();
825  ArrayView<SC> val = valP();
826 
827  LO nStencil = 0; LO tStencil = 0;
828  LO firstCoarseNodeIndex;
829  LO indices[4][3];
830  ia[0] = 0;
831 
832  // Declaration and assignment of fineCoords which holds the coordinates of the fine nodes in 3D.
833  // To do so we pull the nD coordinates from fineCoords and pad the rest with zero vectors...
834  // We also create an nD array that will store the coordinates of the coarse grid nodes.
835  // That array will eventually be used during the construction of coarseCoords, the MultiVector of
836  // coarse nodes coordinates that we store on the coarse level.
837  ArrayRCP< ArrayRCP<double> > lFineCoords(3);
838  ArrayRCP< ArrayRCP<double> > lGhostCoords(3);
839  RCP<const Map> ghostCoordsMap = ghostCoords->getMap();
840  RCP<Xpetra::Vector<double, LO, GO, NO> > zeros = Xpetra::VectorFactory<double, LO, GO, NO>::Build(fineCoordsMap, true);
841  RCP<Xpetra::Vector<double, LO, GO, NO> > ghostZeros = Xpetra::VectorFactory<double, LO, GO, NO>::Build(ghostCoordsMap, true);
842 
843  // Build the MultiVector holding the coarse grid points coordinates.
844  coarseCoords = Xpetra::MultiVectorFactory<double,LO,GO,Node>::Build(coarseCoordsMap, Teuchos::as<size_t>(ndm));
845  ArrayRCP<double> xCoarseNodes; ArrayRCP<double> yCoarseNodes; ArrayRCP<double> zCoarseNodes;
846 
847  if(ndm==1) {
848  lFineCoords[0] = fineCoords->getDataNonConst(0);
849  lFineCoords[1] = zeros->getDataNonConst(0);
850  lFineCoords[2] = zeros->getDataNonConst(0);
851  lGhostCoords[0] = ghostCoords->getDataNonConst(0);
852  lGhostCoords[1] = ghostZeros->getDataNonConst(0);
853  lGhostCoords[2] = ghostZeros->getDataNonConst(0);
854 
855  xCoarseNodes = coarseCoords->getDataNonConst(0);
856  } else if(ndm==2) {
857  lFineCoords[0] = fineCoords->getDataNonConst(0);
858  lFineCoords[1] = fineCoords->getDataNonConst(1);
859  lFineCoords[2] = zeros->getDataNonConst(0);
860  lGhostCoords[0] = ghostCoords->getDataNonConst(0);
861  lGhostCoords[1] = ghostCoords->getDataNonConst(1);
862  lGhostCoords[2] = ghostZeros->getDataNonConst(0);
863 
864  xCoarseNodes= coarseCoords->getDataNonConst(0);
865  yCoarseNodes= coarseCoords->getDataNonConst(1);
866  } else if(ndm==3) {
867  lFineCoords[0] = fineCoords->getDataNonConst(0);
868  lFineCoords[1] = fineCoords->getDataNonConst(1);
869  lFineCoords[2] = fineCoords->getDataNonConst(2);
870  lGhostCoords[0] = ghostCoords->getDataNonConst(0);
871  lGhostCoords[1] = ghostCoords->getDataNonConst(1);
872  lGhostCoords[2] = ghostCoords->getDataNonConst(2);
873 
874  xCoarseNodes = coarseCoords->getDataNonConst(0);
875  yCoarseNodes = coarseCoords->getDataNonConst(1);
876  zCoarseNodes = coarseCoords->getDataNonConst(2);
877  }
878 
879  // Finally let's load the PIDs of the coarse nodes on the coarse grid
880  // to later perform permutation of stencil entries to conform with
881  // Tpetra matrix ordering when calling setAllValues.
882 
883  // Loop over the fine nodes and compute their interpolation stencils
884  // Get the rank of the local process to order entries by PIDs in the
885  // prolongator matrix rows
886  GO currentCoarseNode = 0;
887  for(LO i = 0; i < lNumFineNodes; ++i) {
888 
889  // Get point indices on fine grid
890  {
891  std::div_t tmp;
892  tmp=std::div(i,lFineNodesPerDir[0]*lFineNodesPerDir[1]);
893  indices[0][2] = tmp.quot;
894  tmp=std::div(tmp.rem,lFineNodesPerDir[0]);
895  indices[0][1] = tmp.quot;
896  indices[0][0] = tmp.rem;
897 
898  // Get ref point indices on coarse grid
899  tmp=std::div(indices[0][0] - offsets[0],coarseRate[0]);
900  indices[1][0] = tmp.quot;
901  tmp=std::div(indices[0][1] - offsets[1],coarseRate[1]);
902  indices[1][1] = tmp.quot;
903  tmp=std::div(indices[0][2] - offsets[2],coarseRate[2]);
904  indices[1][2] = tmp.quot;
905  }
906 
907  // location "flags" indicate if the current node is on a coarse
908  // face, edge or node.
909  indices[2][0] = (indices[0][0] + offsets[0]) % coarseRate[0];
910  indices[2][1] = (indices[0][1] + offsets[1]) % coarseRate[1];
911  indices[2][2] = (indices[0][2] + offsets[2]) % coarseRate[2];
912 
913  // Get indices of ref point on fine grid
914  indices[3][0] = indices[1][0]*coarseRate[0];
915  indices[3][1] = indices[1][1]*coarseRate[1];
916  indices[3][2] = indices[1][2]*coarseRate[2];
917 
918  if( (indices[0][0] == lFineNodesPerDir[0]-1) && !ghostInterface[1] ) {
919  indices[1][0] = lCoarseNodesPerDir[0]-1;
920  indices[2][0] = 0;
921  indices[3][0] = lFineNodesPerDir[0]-1;
922  }
923  if( (indices[0][1] == lFineNodesPerDir[1]-1) && !ghostInterface[3] ) {
924  indices[1][1] = lCoarseNodesPerDir[1]-1;
925  indices[2][1] = 0;
926  indices[3][1] = lFineNodesPerDir[1]-1;
927  }
928  if( (indices[0][2] == lFineNodesPerDir[2]-1) && !ghostInterface[5] ) {
929  indices[1][2] = lCoarseNodesPerDir[2]-1;
930  indices[2][2] = 0;
931  indices[3][2] = lFineNodesPerDir[2]-1;
932  }
933 
934  Array<GO> currentNodeIndices(3);
935  currentNodeIndices[0] = gStartIndices[0] + Teuchos::as<GO>(indices[0][mapDirL2G[0]]);
936  currentNodeIndices[1] = gStartIndices[1] + Teuchos::as<GO>(indices[0][mapDirL2G[1]]);
937  currentNodeIndices[2] = gStartIndices[2] + Teuchos::as<GO>(indices[0][mapDirL2G[2]]);
938 
939  // Assuming lexicographic indexing the coarse nodes forming a prism
940  // around fine node "i" are selected and store them in connec.
941  // Two tricky things to be careful about:
942  // - are we using coarseRate or endRate?
943  // --> check indices and set rate correctly
944  // - are we on the east, north or top face?
945  // --> if so fix firstCoarseNodeIndex to make sure there is
946  // a node to the east, north and top of firstCoarseNodeIndex
947  LO rate[3];
948  if(currentNodeIndices[0] >= gFineNodesPerDir[0] - endRate[0] - 1) {
949  rate[0] = endRate[0];
950  } else {
951  rate[0] = coarseRate[0];
952  }
953  if(currentNodeIndices[1] >= gFineNodesPerDir[1] - endRate[1] - 1) {
954  rate[1] = endRate[1];
955  } else {
956  rate[1] = coarseRate[1];
957  }
958  if(currentNodeIndices[2] >= gFineNodesPerDir[2] - endRate[2] - 1) {
959  rate[2] = endRate[2];
960  } else {
961  rate[2] = coarseRate[2];
962  }
963  if(ndm < 3) { rate[2] = 0;}
964  if(ndm < 2) { rate[1] = 0;}
965 
966  // We need to check whether we are on the edge of the mesh in which case we need to adjust the coarse nodes
967  firstCoarseNodeIndex = 0;
968  Array<GO> firstCoarseNodeIndices(3); // These are fine grid indices, divide by coarseRate[i] to get coarse grid indices
969  if((currentNodeIndices[2] == gFineNodesPerDir[2] -1) && (endRate[2] == coarseRate[2])) {
970  // If we are on the last node and have a endRate == coarseRate we need to take the coarse node below the current node
971  firstCoarseNodeIndices[2] = ((currentNodeIndices[2] / coarseRate[2]) - 1) * coarseRate[2];
972  } else {
973  firstCoarseNodeIndices[2] = (currentNodeIndices[2] / coarseRate[2]) * coarseRate[2];
974  }
975  if((currentNodeIndices[1] == gFineNodesPerDir[1] -1) && (endRate[1] == coarseRate[1])) {
976  firstCoarseNodeIndices[1] = ((currentNodeIndices[1] / coarseRate[1]) - 1) * coarseRate[1];
977  } else {
978  firstCoarseNodeIndices[1] = (currentNodeIndices[1] / coarseRate[1]) * coarseRate[1];
979  }
980  if((currentNodeIndices[0] == gFineNodesPerDir[0] -1) && (endRate[0] == coarseRate[0])) {
981  firstCoarseNodeIndices[0] = ((currentNodeIndices[0] / coarseRate[0]) - 1) * coarseRate[0];
982  } else {
983  firstCoarseNodeIndices[0] = (currentNodeIndices[0] / coarseRate[0]) * coarseRate[0];
984  }
985  firstCoarseNodeIndex += firstCoarseNodeIndices[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0]
986  + firstCoarseNodeIndices[1]*gFineNodesPerDir[0] + firstCoarseNodeIndices[0];
987 
988  GO firstCoarseNodeOnCoarseGridIndex;
989  {
990  GO tmpInds[4];
991  tmpInds[2] = firstCoarseNodeIndex / fineNodesPerCoarseSlab;
992  tmpInds[3] = firstCoarseNodeIndex % fineNodesPerCoarseSlab;
993  tmpInds[1] = tmpInds[3] / fineNodesPerCoarsePlane;
994  tmpInds[0] = (tmpInds[3] % fineNodesPerCoarsePlane) / coarseRate[0];
995  firstCoarseNodeOnCoarseGridIndex = tmpInds[2]*coarseNodesPerCoarseLayer + tmpInds[1]*gCoarseNodesPerDir[0] + tmpInds[0];
996  }
997 
998  LO coarseGhosts[8] = {0, 0, 0, 0, 0, 0, 0, 0};
999  if( ghostInterface[0] ) {
1000  if( ((indices[0][mapDirG2L[0]] < rate[0] - offsets[0]) && offsets[0] != 0)
1001  || (((currentNodeIndices[0] == gFineNodesPerDir[0] -1) && (indices[0][mapDirG2L[0]] < rate[0] - offsets[0] + 1)) && offsets[0] != 0)
1002  || ((currentNodeIndices[0] == gFineNodesPerDir[0] -1) && lFineNodesPerDir[mapDirG2L[0]] == 1) ) {
1003  coarseGhosts[0] = 1;
1004  coarseGhosts[2] = 1;
1005  coarseGhosts[4] = 1;
1006  coarseGhosts[6] = 1;
1007  }
1008  }
1009  if( ghostInterface[1] && (indices[0][mapDirG2L[0]] > lFineNodesPerDir[mapDirG2L[0]] - offsets[3] - 2) ) {
1010  coarseGhosts[1] = 1;
1011  coarseGhosts[3] = 1;
1012  coarseGhosts[5] = 1;
1013  coarseGhosts[7] = 1;
1014  }
1015  if( ghostInterface[2] ) {
1016  if( ((indices[0][mapDirG2L[1]] < rate[1] - offsets[1]) && offsets[1] != 0)
1017  || (((currentNodeIndices[1] == gFineNodesPerDir[1] -1) && (indices[0][mapDirG2L[1]] < rate[1] - offsets[1] + 1)) && offsets[1] != 0)
1018  || ((currentNodeIndices[1] == gFineNodesPerDir[1] -1) && lFineNodesPerDir[mapDirG2L[1]] == 1) ) {
1019  coarseGhosts[0] = 1;
1020  coarseGhosts[1] = 1;
1021  coarseGhosts[4] = 1;
1022  coarseGhosts[5] = 1;
1023  }
1024  }
1025  if( ghostInterface[3] && (indices[0][mapDirG2L[1]] > lFineNodesPerDir[mapDirG2L[1]] - offsets[4] - 2) ) {
1026  coarseGhosts[2] = 1;
1027  coarseGhosts[3] = 1;
1028  coarseGhosts[6] = 1;
1029  coarseGhosts[7] = 1;
1030  }
1031  if( ghostInterface[4] ) {
1032  if( ((indices[0][mapDirG2L[2]] < rate[2] - offsets[2]) && offsets[2] != 0)
1033  || (((currentNodeIndices[2] == gFineNodesPerDir[2] -1) && (indices[0][mapDirG2L[2]] < rate[2] - offsets[2] + 1)) && offsets[2] != 0)
1034  || ((currentNodeIndices[2] == gFineNodesPerDir[2] -1) && lFineNodesPerDir[mapDirG2L[2]] == 1) ) {
1035  coarseGhosts[0] = 1;
1036  coarseGhosts[1] = 1;
1037  coarseGhosts[2] = 1;
1038  coarseGhosts[3] = 1;
1039  }
1040  }
1041  if( ghostInterface[5] && (indices[0][mapDirG2L[2]] > lFineNodesPerDir[mapDirG2L[2]] - offsets[5] - 2) ) {
1042  coarseGhosts[4] = 1;
1043  coarseGhosts[5] = 1;
1044  coarseGhosts[6] = 1;
1045  coarseGhosts[7] = 1;
1046  }
1047 
1048  GO firstGhostNodeIndices[3], firstGhostNodeIndex;
1049  if(currentNodeIndices[0] == gFineNodesPerDir[0] - 1) {
1050  firstGhostNodeIndices[0] = (currentNodeIndices[0]-rate[0]) - (currentNodeIndices[0]-rate[0])%coarseRate[0];
1051  } else {
1052  firstGhostNodeIndices[0] = currentNodeIndices[0] - currentNodeIndices[0]%coarseRate[0];
1053  }
1054  if(currentNodeIndices[1] == gFineNodesPerDir[1] - 1) {
1055  firstGhostNodeIndices[1] = (currentNodeIndices[1]-rate[1]) - (currentNodeIndices[1]-rate[1])%coarseRate[1];
1056  } else {
1057  firstGhostNodeIndices[1] = currentNodeIndices[1] - currentNodeIndices[1]%coarseRate[1];
1058  }
1059  if(currentNodeIndices[2] == gFineNodesPerDir[2] - 1) {
1060  firstGhostNodeIndices[2] = (currentNodeIndices[2]-rate[2]) - (currentNodeIndices[2]-rate[2])%coarseRate[2];
1061  } else {
1062  firstGhostNodeIndices[2] = currentNodeIndices[2] - currentNodeIndices[2]%coarseRate[2];
1063  }
1064  firstGhostNodeIndex = firstGhostNodeIndices[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0] + firstGhostNodeIndices[1]*gFineNodesPerDir[0] + firstGhostNodeIndices[0];
1065 
1066  Array<GO> gCoarseNodesOnCoarseGridIndices(8);
1067  GO ghostNodeIndex, ghostNodeOnCoarseGridIndex;
1068  GO coarseNodeIndex, coarseNodeOnCoarseGridIndex;
1069  LO ind;
1070  double connec[9][3];
1071  LO currentGhostLID;
1072  Array<LO> connecPIDs(8); // A PID of -1 indicate a local element.
1073  Array<GO> connecLGIDs(8);
1074  // First load the coordinates of the node of interest
1075  for(LO dim = 0; dim < 3; ++dim) {connec[0][dim] = lFineCoords[dim][i];}
1076  // Then loop over the three spatial dimension and load the 8 coarse nodes
1077  // required for the linear interpolation.
1078  for(LO ind2 = 0; ind2 < 2; ++ind2) {
1079  for(LO ind1 = 0; ind1 < 2; ++ind1) {
1080  for(LO ind0 = 0; ind0 < 2; ++ind0) {
1081  ind = 4*ind2 + 2*ind1 + ind0;
1082  // Check whether ghost nodes are needed for the current fine point.
1083  if(coarseGhosts[ind] == 1) {
1084  // Get the global ghost node index and load its coordinates
1085  ghostNodeIndex = firstGhostNodeIndex + Teuchos::as<GO>(ind2*rate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0]
1086  + ind1*rate[1]*gFineNodesPerDir[0]
1087  + ind0*rate[0]);
1088  ghostNodeOnCoarseGridIndex = firstCoarseNodeOnCoarseGridIndex + Teuchos::as<GO>(ind2*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0]
1089  + ind1*gCoarseNodesPerDir[0]
1090  + ind0);
1091  currentGhostLID = ghostMap->getLocalElement(ghostNodeIndex);
1092  connecPIDs[ind] = ghostsPIDs[currentGhostLID];
1093  connecLGIDs[ind] = ghostNodeIndex;
1094  for(LO dim = 0; dim < 3; ++dim) {connec[ind + 1][dim] = lGhostCoords[dim][currentGhostLID];}
1095  gCoarseNodesOnCoarseGridIndices[4*ind2 + 2*ind1 + ind0] = ghostNodeOnCoarseGridIndex;
1096  } else {
1097  // Get the local coarse node index and load its coordinates
1098  coarseNodeIndex = firstCoarseNodeIndex + Teuchos::as<GO>(ind2*rate[2]*gFineNodesPerDir[1]*gFineNodesPerDir[0]
1099  + ind1*rate[1]*gFineNodesPerDir[0]
1100  + ind0*rate[0]);
1101  coarseNodeOnCoarseGridIndex = firstCoarseNodeOnCoarseGridIndex + Teuchos::as<GO>(ind2*gCoarseNodesPerDir[1]*gCoarseNodesPerDir[0]
1102  + ind1*gCoarseNodesPerDir[0]
1103  + ind0);
1104  for(LO dim = 0; dim < 3; ++dim) {connec[ind + 1][dim] = lFineCoords[dim][fineCoordsMap->getLocalElement(coarseNodeIndex)];}
1105  gCoarseNodesOnCoarseGridIndices[4*ind2 + 2*ind1 + ind0] = coarseNodeOnCoarseGridIndex;
1106  connecPIDs[ind] = -1;
1107  connecLGIDs[ind] = coarseCoordsMap->getLocalElement(coarseNodeOnCoarseGridIndex);
1108  }
1109  }
1110  }
1111  }
1112 
1113  // Compute the actual geometric interpolation stencil
1114  SC stencil[8] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
1115  ComputeStencil(ndm, currentNodeIndices, firstCoarseNodeIndices, rate, connec, interpolationOrder, stencil);
1116 
1117  // Finally check whether the fine node is on a coarse: node, edge or face
1118  // and select accordingly the non-zero values from the stencil and the
1119  // corresponding column indices
1120  LO nzIndStencil[8] = {0,0,0,0,0,0,0,0};
1121  if( ((currentNodeIndices[0] % coarseRate[0] == 0) || currentNodeIndices[0] == gFineNodesPerDir[0] - 1)
1122  && ((currentNodeIndices[1] % coarseRate[1] == 0) || currentNodeIndices[1] == gFineNodesPerDir[1] - 1)
1123  && ((currentNodeIndices[2] % coarseRate[2] == 0) || currentNodeIndices[2] == gFineNodesPerDir[2] - 1) ) {
1124  if(ndm==1) {
1125  xCoarseNodes[currentCoarseNode] = lFineCoords[0][i];
1126  } else if(ndm==2) {
1127  xCoarseNodes[currentCoarseNode] = lFineCoords[0][i];
1128  yCoarseNodes[currentCoarseNode] = lFineCoords[1][i];
1129  } else if(ndm==3) {
1130  xCoarseNodes[currentCoarseNode] = lFineCoords[0][i];
1131  yCoarseNodes[currentCoarseNode] = lFineCoords[1][i];
1132  zCoarseNodes[currentCoarseNode] = lFineCoords[2][i];
1133  }
1134  ++currentCoarseNode;
1135 
1136  if(currentNodeIndices[0] == gFineNodesPerDir[0] - 1) {
1137  nzIndStencil[0] += 1;
1138  }
1139  if((currentNodeIndices[1] == gFineNodesPerDir[1] - 1) && (ndm > 1)) {
1140  nzIndStencil[0] += 2;
1141  }
1142  if((currentNodeIndices[2] == gFineNodesPerDir[2] - 1) && (ndm > 2)) {
1143  nzIndStencil[0] += 4;
1144  }
1145  nStencil = 1;
1146  for(LO k = 0; k < 8; ++k) {
1147  nzIndStencil[k] = nzIndStencil[0];
1148  }
1149  } else {
1150  nStencil = 8;
1151  for(LO k = 0; k < 8; ++k) {
1152  nzIndStencil[k] = k;
1153  }
1154  }
1155 
1156  // Here the values are filled in the Crs matrix arrays
1157  // This is basically the only place these variables are modified
1158  // Hopefully this makes handling system of PDEs easy!
1159 
1160  // Loop on dofsPerNode and process each row for the current Node
1161 
1162  Array<LO> permutation(8);
1163  for(LO k = 0; k < 8; ++k) {permutation[k] = k;}
1164 
1165  // Sort nodes by PIDs using stable sort to keep sublist ordered by LIDs and GIDs
1166  sh_sort2(connecPIDs.begin(),connecPIDs.end(), permutation.begin(), permutation.end());
1167 
1168  for(LO j = 0; j < dofsPerNode; ++j) {
1169  ia[i*dofsPerNode + j + 1] = ia[i*dofsPerNode + j] + nStencil;
1170  for(LO k = 0; k < nStencil; ++k) {
1171  ja [ia[i*dofsPerNode + j] + k] = colMapP->getLocalElement(gCoarseNodesOnCoarseGridIndices[nzIndStencil[permutation[k]]]*dofsPerNode + j);
1172  val[ia[i*dofsPerNode + j] + k] = stencil[nzIndStencil[permutation[k]]];
1173  }
1174  // Add the stencil for each degree of freedom.
1175  tStencil += nStencil;
1176  }
1177  } // End loop over fine nodes
1178 
1179  if (rowMapP->lib() == Xpetra::UseTpetra) {
1180  // - Cannot resize for Epetra, as it checks for same pointers
1181  // - Need to resize for Tpetra, as it check ().size() == ia[numRows]
1182  // NOTE: these invalidate ja and val views
1183  jaP .resize(tStencil);
1184  valP.resize(tStencil);
1185  }
1186 
1187  // Set the values of the prolongation operators into the CrsMatrix P and call FillComplete
1188  PCrs->setAllValues(iaP, jaP, valP);
1189  PCrs->expertStaticFillComplete(domainMapP,rowMapP);
1190 
1191  } // MakeGeneralGeometricP
1192 
1193  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1194  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeStencil(
1195  const LO numDimensions, const Array<GO> currentNodeIndices,
1196  const Array<GO> coarseNodeIndices, const LO rate[3],
1197  const double coord[9][3], const int interpolationOrder,
1198  SC stencil[8]) const {
1199 
1200  TEUCHOS_TEST_FOR_EXCEPTION((interpolationOrder > 1) || (interpolationOrder < 0),
1201  Exceptions::RuntimeError,
1202  "The interpolation order can be set to 0 or 1 only.");
1203 
1204  if(interpolationOrder == 0) {
1205  ComputeConstantInterpolationStencil(numDimensions, currentNodeIndices, coarseNodeIndices, rate, stencil);
1206  } else if(interpolationOrder == 1) {
1207  ComputeLinearInterpolationStencil(numDimensions, coord, stencil);
1208  }
1209 
1210  } // End ComputeStencil
1211 
1212  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1213  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeConstantInterpolationStencil(
1214  const LO numDimensions, const Array<GO> currentNodeIndices,
1215  const Array<GO> coarseNodeIndices, const LO rate[3], SC stencil[8]) const {
1216 
1217  LO coarseNode = 0;
1218  if(numDimensions > 2) {
1219  if((currentNodeIndices[2] - coarseNodeIndices[2]) > (rate[2] / 2)) {
1220  coarseNode += 4;
1221  }
1222  }
1223  if(numDimensions > 1) {
1224  if((currentNodeIndices[1] - coarseNodeIndices[1]) > (rate[1] / 2)) {
1225  coarseNode += 2;
1226  }
1227  }
1228  if((currentNodeIndices[0] - coarseNodeIndices[0]) > (rate[0] / 2)) {
1229  coarseNode += 1;
1230  }
1231  stencil[coarseNode] = 1.0;
1232 
1233  }
1234 
1235  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1236  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::ComputeLinearInterpolationStencil(
1237  const LO numDimensions, const double coord[9][3], SC stencil[8]) const {
1238 
1239  // 7 8 Find xi, eta and zeta such that
1240  // x---------x
1241  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
1242  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
1243  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
1244  // | | *P | |
1245  // | x------|--x We can do this with a Newton solver:
1246  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
1247  // |/ |/ We compute the Jacobian and iterate until convergence...
1248  // z y x---------x
1249  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
1250  // |/ give us the weights for the interpolation stencil!
1251  // o---x
1252  //
1253 
1254  Teuchos::SerialDenseMatrix<LO,double> Jacobian(numDimensions, numDimensions);
1255  Teuchos::SerialDenseVector<LO,double> residual(numDimensions);
1256  Teuchos::SerialDenseVector<LO,double> solutionDirection(numDimensions);
1257  Teuchos::SerialDenseVector<LO,double> paramCoords(numDimensions);
1258  Teuchos::SerialDenseSolver<LO,double> problem;
1259  LO numTerms = std::pow(2,numDimensions), iter = 0, max_iter = 5;
1260  double functions[4][8], norm_ref = 1, norm2 = 1, tol = 1e-5;
1261  paramCoords.size(numDimensions);
1262 
1263  while( (iter < max_iter) && (norm2 > tol*norm_ref) ) {
1264  ++iter;
1265  norm2 = 0;
1266  solutionDirection.size(numDimensions);
1267  residual.size(numDimensions);
1268  Jacobian = 0.0;
1269 
1270  // Compute Jacobian and Residual
1271  GetInterpolationFunctions(numDimensions, paramCoords, functions);
1272  for(LO i = 0; i < numDimensions; ++i) {
1273  residual(i) = coord[0][i]; // Add coordinates from point of interest
1274  for(LO k = 0; k < numTerms; ++k) {
1275  residual(i) -= functions[0][k]*coord[k+1][i]; // Remove contribution from all coarse points
1276  }
1277  if(iter == 1) {
1278  norm_ref += residual(i)*residual(i);
1279  if(i == numDimensions - 1) {
1280  norm_ref = std::sqrt(norm_ref);
1281  }
1282  }
1283 
1284  for(LO j = 0; j < numDimensions; ++j) {
1285  for(LO k = 0; k < numTerms; ++k) {
1286  Jacobian(i,j) += functions[j+1][k]*coord[k+1][i];
1287  }
1288  }
1289  }
1290 
1291  // Set Jacobian, Vectors and solve problem
1292  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
1293  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
1294  problem.factorWithEquilibration(true);
1295  problem.solve();
1296  problem.unequilibrateLHS();
1297 
1298  for(LO i = 0; i < numDimensions; ++i) {
1299  paramCoords(i) = paramCoords(i) + solutionDirection(i);
1300  }
1301 
1302  // Recompute Residual norm
1303  GetInterpolationFunctions(numDimensions, paramCoords, functions); // These need to be recomputed with new paramCoords!
1304  for(LO i = 0; i < numDimensions; ++i) {
1305  double tmp = coord[0][i];
1306  for(LO k = 0; k < numTerms; ++k) {
1307  tmp -= functions[0][k]*coord[k+1][i];
1308  }
1309  norm2 += tmp*tmp;
1310  tmp = 0;
1311  }
1312  norm2 = std::sqrt(norm2);
1313  }
1314 
1315  // Load the interpolation values onto the stencil.
1316  for(LO i = 0; i < 8; ++i) {
1317  stencil[i] = functions[0][i];
1318  }
1319 
1320  }
1321 
1322  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1323  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector<LO,double> parameters, double functions[4][8]) const {
1324  double xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
1325  if(numDimensions == 1) {
1326  xi = parameters[0];
1327  denominator = 2.0;
1328  } else if(numDimensions == 2) {
1329  xi = parameters[0];
1330  eta = parameters[1];
1331  denominator = 4.0;
1332  } else if(numDimensions == 3) {
1333  xi = parameters[0];
1334  eta = parameters[1];
1335  zeta = parameters[2];
1336  denominator = 8.0;
1337  }
1338 
1339  functions[0][0] = (1.0 - xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
1340  functions[0][1] = (1.0 + xi)*(1.0 - eta)*(1.0 - zeta) / denominator;
1341  functions[0][2] = (1.0 - xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
1342  functions[0][3] = (1.0 + xi)*(1.0 + eta)*(1.0 - zeta) / denominator;
1343  functions[0][4] = (1.0 - xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
1344  functions[0][5] = (1.0 + xi)*(1.0 - eta)*(1.0 + zeta) / denominator;
1345  functions[0][6] = (1.0 - xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
1346  functions[0][7] = (1.0 + xi)*(1.0 + eta)*(1.0 + zeta) / denominator;
1347 
1348  functions[1][0] = -(1.0 - eta)*(1.0 - zeta) / denominator;
1349  functions[1][1] = (1.0 - eta)*(1.0 - zeta) / denominator;
1350  functions[1][2] = -(1.0 + eta)*(1.0 - zeta) / denominator;
1351  functions[1][3] = (1.0 + eta)*(1.0 - zeta) / denominator;
1352  functions[1][4] = -(1.0 - eta)*(1.0 + zeta) / denominator;
1353  functions[1][5] = (1.0 - eta)*(1.0 + zeta) / denominator;
1354  functions[1][6] = -(1.0 + eta)*(1.0 + zeta) / denominator;
1355  functions[1][7] = (1.0 + eta)*(1.0 + zeta) / denominator;
1356 
1357  functions[2][0] = -(1.0 - xi)*(1.0 - zeta) / denominator;
1358  functions[2][1] = -(1.0 + xi)*(1.0 - zeta) / denominator;
1359  functions[2][2] = (1.0 - xi)*(1.0 - zeta) / denominator;
1360  functions[2][3] = (1.0 + xi)*(1.0 - zeta) / denominator;
1361  functions[2][4] = -(1.0 - xi)*(1.0 + zeta) / denominator;
1362  functions[2][5] = -(1.0 + xi)*(1.0 + zeta) / denominator;
1363  functions[2][6] = (1.0 - xi)*(1.0 + zeta) / denominator;
1364  functions[2][7] = (1.0 + xi)*(1.0 + zeta) / denominator;
1365 
1366  functions[3][0] = -(1.0 - xi)*(1.0 - eta) / denominator;
1367  functions[3][1] = -(1.0 + xi)*(1.0 - eta) / denominator;
1368  functions[3][2] = -(1.0 - xi)*(1.0 + eta) / denominator;
1369  functions[3][3] = -(1.0 + xi)*(1.0 + eta) / denominator;
1370  functions[3][4] = (1.0 - xi)*(1.0 - eta) / denominator;
1371  functions[3][5] = (1.0 + xi)*(1.0 - eta) / denominator;
1372  functions[3][6] = (1.0 - xi)*(1.0 + eta) / denominator;
1373  functions[3][7] = (1.0 + xi)*(1.0 + eta) / denominator;
1374 
1375  } // End GetInterpolationFunctions
1376 
1377  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1378  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort_permute(
1379  const typename Teuchos::Array<LocalOrdinal>::iterator& first1,
1380  const typename Teuchos::Array<LocalOrdinal>::iterator& last1,
1381  const typename Teuchos::Array<LocalOrdinal>::iterator& first2,
1382  const typename Teuchos::Array<LocalOrdinal>::iterator& last2) const
1383  {
1384  typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1385  DT n = last1 - first1;
1386  DT m = n / 2;
1387  DT z = Teuchos::OrdinalTraits<DT>::zero();
1388  while (m > z)
1389  {
1390  DT max = n - m;
1391  for (DT j = 0; j < max; j++)
1392  {
1393  for (DT k = j; k >= 0; k-=m)
1394  {
1395  if (first1[first2[k+m]] >= first1[first2[k]])
1396  break;
1397  std::swap(first2[k+m], first2[k]);
1398  }
1399  }
1400  m = m/2;
1401  }
1402  }
1403 
1404  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
1405  void GeneralGeometricPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::sh_sort2(
1406  const typename Teuchos::Array<LocalOrdinal>::iterator& first1,
1407  const typename Teuchos::Array<LocalOrdinal>::iterator& last1,
1408  const typename Teuchos::Array<LocalOrdinal>::iterator& first2,
1409  const typename Teuchos::Array<LocalOrdinal>::iterator& last2) const
1410  {
1411  typedef typename std::iterator_traits<typename Teuchos::Array<LocalOrdinal>::iterator>::difference_type DT;
1412  DT n = last1 - first1;
1413  DT m = n / 2;
1414  DT z = Teuchos::OrdinalTraits<DT>::zero();
1415  while (m > z)
1416  {
1417  DT max = n - m;
1418  for (DT j = 0; j < max; j++)
1419  {
1420  for (DT k = j; k >= 0; k-=m)
1421  {
1422  if (first1[k+m] >= first1[k])
1423  break;
1424  std::swap(first1[k+m], first1[k]);
1425  std::swap(first2[k+m], first2[k]);
1426  }
1427  }
1428  m = m/2;
1429  }
1430  }
1431 
1432 } //namespace MueLu
1433 
1434 #define MUELU_GENERALGEOMETRICPFACTORY_SHORT
1435 #endif // MUELU_GENERALGEOMETRICPFACTORY_DEF_HPP
MueLu::GeneralGeometricPFactory::sh_sort2
void sh_sort2(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1405
MueLu_FactoryManagerBase.hpp
MueLu::GeneralGeometricPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:76
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:191
MueLu::GeneralGeometricPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const LO numDimension, const double coord[9][3], SC stencil[8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1236
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:76
MueLu::GeneralGeometricPFactory::sh_sort_permute
void sh_sort_permute(const typename Teuchos::Array< LocalOrdinal >::iterator &first1, const typename Teuchos::Array< LocalOrdinal >::iterator &last1, const typename Teuchos::Array< LocalOrdinal >::iterator &first2, const typename Teuchos::Array< LocalOrdinal >::iterator &last2) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1378
MueLu::GeneralGeometricPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:132
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu_MasterList.hpp
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
MueLu::GeneralGeometricPFactory::ComputeStencil
void ComputeStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], const double coord[9][3], const int interpolationOrder, SC stencil[8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1194
MueLu::GeneralGeometricPFactory::MakeGeneralGeometricP
void MakeGeneralGeometricP(LO const numDimension, const Array< LO > mapDirL2G, const Array< LO > mapDirG2L, const Array< LO > lFineNodesPerDir, const Array< LO > lCoarseNodesPerDir, Array< GO > gCoarseNodesPerDir, Array< GO > gFineNodesPerDir, ArrayRCP< LO > const coarseRate, LO const endRate[3], LO const offsets[6], bool const ghostInterface[6], const RCP< Xpetra::MultiVector< double, LO, GO, NO > > &fCoords, LO const nnzP, LO const dofsPerNode, RCP< const Map > &stridedDomainMapP, RCP< Matrix > &Amat, RCP< Matrix > &P, RCP< Xpetra::MultiVector< double, LO, GO, NO > > &cCoords, Array< GO > ghostsGIDs, int interpolationOrder) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:587
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::GeneralGeometricPFactory::ComputeConstantInterpolationStencil
void ComputeConstantInterpolationStencil(const LO numDimension, const Array< GO > currentNodeIndices, const Array< GO > coarseNodeIndices, const LO rate[3], SC stencil[8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1213
MueLu_GeneralGeometricPFactory_decl.hpp
MueLu::GeneralGeometricPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimension, const Teuchos::SerialDenseVector< LO, double > parameters, double functions[4][8]) const
Definition: MueLu_GeneralGeometricPFactory_def.hpp:1323
MueLu_Monitor.hpp
MueLu::Level::IsAvailable
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need&#39;s value has been saved.
Definition: MueLu_Level.hpp:332
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:70
MueLu::GeneralGeometricPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:127
MueLu::GeneralGeometricPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeneralGeometricPFactory_def.hpp:95
MueLu::Level::DeclareInput
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
Definition: MueLu_Level.cpp:150