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