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