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