46 #ifndef MUELU_REGIONRFACTORY_DEF_HPP 47 #define MUELU_REGIONRFACTORY_DEF_HPP 49 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 51 #include <Xpetra_Matrix.hpp> 52 #include <Xpetra_CrsGraphFactory.hpp> 63 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
64 RCP<const ParameterList> RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
65 GetValidParameterList()
const {
66 RCP<ParameterList> validParamList = rcp(
new ParameterList());
68 validParamList->set<RCP<const FactoryBase> >(
"A", Teuchos::null,
69 "Generating factory of the matrix A");
70 validParamList->set<RCP<const FactoryBase> >(
"numDimensions", Teuchos::null,
71 "Number of spatial dimensions in the problem.");
72 validParamList->set<RCP<const FactoryBase> >(
"lNodesPerDim", Teuchos::null,
73 "Number of local nodes per spatial dimension on the fine grid.");
74 validParamList->set<RCP<const FactoryBase> >(
"Nullspace", Teuchos::null,
75 "Fine level nullspace used to construct the coarse level nullspace.");
76 validParamList->set<RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
77 "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
78 validParamList->set<
bool> (
"keep coarse coords",
false,
"Flag to keep coordinates for special coarse grid solve");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
85 DeclareInput(Level& fineLevel, Level& )
const {
87 Input(fineLevel,
"A");
88 Input(fineLevel,
"numDimensions");
89 Input(fineLevel,
"lNodesPerDim");
90 Input(fineLevel,
"Nullspace");
91 Input(fineLevel,
"Coordinates");
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
97 Build(Level& fineLevel, Level& coarseLevel)
const {
100 RCP<Teuchos::FancyOStream> out;
101 if(
const char* dbg = std::getenv(
"MUELU_REGIONRFACTORY_DEBUG")) {
102 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
103 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
105 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
108 *out <<
"Starting RegionRFactory::Build." << std::endl;
111 const int numDimensions = Get<int>(fineLevel,
"numDimensions");
112 Array<LO> lFineNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
114 Array<LO> lNodesPerDim = Get<Array<LO> >(fineLevel,
"lNodesPerDim");
115 for(
int dim = 0; dim < numDimensions; ++dim) {
116 lFineNodesPerDim[dim] = lNodesPerDim[dim];
119 *out <<
"numDimensions " << numDimensions <<
" and lFineNodesPerDim: " << lFineNodesPerDim
123 if(numDimensions < 1 || numDimensions > 3) {
124 throw std::runtime_error(
"numDimensions must be 1, 2 or 3!");
126 for(
int dim = 0; dim < numDimensions; ++dim) {
127 if(lFineNodesPerDim[dim] % 3 != 1) {
128 throw std::runtime_error(
"The number of fine node in each direction need to be 3n+1");
131 Array<LO> lCoarseNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
133 const RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel,
"A");
135 RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
136 fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel,
"Coordinates");
137 if(static_cast<int>(fineCoordinates->getNumVectors()) != numDimensions) {
138 throw std::runtime_error(
"The number of vectors in the coordinates is not equal to numDimensions!");
146 if(numDimensions == 1) {
147 throw std::runtime_error(
"RegionRFactory no implemented for 1D case yet.");
148 }
else if(numDimensions == 2) {
149 throw std::runtime_error(
"RegionRFactory no implemented for 2D case yet.");
150 }
else if(numDimensions == 3) {
151 Build3D(numDimensions, lFineNodesPerDim, A, fineCoordinates,
152 R, coarseCoordinates, lCoarseNodesPerDim);
155 const Teuchos::ParameterList& pL = GetParameterList();
158 RCP<ParameterList> Tparams;
159 if(pL.isSublist(
"matrixmatrix: kernel params"))
160 Tparams=rcp(
new ParameterList(pL.sublist(
"matrixmatrix: kernel params")));
162 Tparams= rcp(
new ParameterList);
165 *out <<
"Compute P=R^t" << std::endl;
167 Tparams->set(
"compute global constants: temporaries",Tparams->get(
"compute global constants: temporaries",
false));
168 Tparams->set(
"compute global constants", Tparams->get(
"compute global constants",
false));
169 std::string label =
"MueLu::RegionR-transR" +
Teuchos::toString(coarseLevel.GetLevelID());
172 *out <<
"Compute coarse nullspace" << std::endl;
173 RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel,
"Nullspace");
174 RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(R->getRowMap(),
175 fineNullspace->getNumVectors());
176 R->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
177 Teuchos::ScalarTraits<SC>::zero());
179 *out <<
"Set data on coarse level" << std::endl;
180 Set(coarseLevel,
"numDimensions", numDimensions);
181 Set(coarseLevel,
"lNodesPerDim", lCoarseNodesPerDim);
182 Set(coarseLevel,
"Nullspace", coarseNullspace);
183 Set(coarseLevel,
"Coordinates", coarseCoordinates);
184 if(pL.get<
bool>(
"keep coarse coords")) {
185 coarseLevel.Set<RCP<realvaluedmultivector_type> >(
"Coordinates2", coarseCoordinates,
NoFactory::get());
188 R->SetFixedBlockSize(A->GetFixedBlockSize());
189 P->SetFixedBlockSize(A->GetFixedBlockSize());
191 Set(coarseLevel,
"R", R);
192 Set(coarseLevel,
"P", P);
196 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
197 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
198 Build3D(
const int numDimensions,
199 Teuchos::Array<LocalOrdinal>& lFineNodesPerDim,
200 const RCP<Matrix>& A,
204 Teuchos::Array<LocalOrdinal>& lCoarseNodesPerDim)
const {
205 using local_matrix_type =
typename CrsMatrix::local_matrix_type;
206 using local_graph_type =
typename local_matrix_type::staticcrsgraph_type;
207 using row_map_type =
typename local_matrix_type::row_map_type::non_const_type;
208 using entries_type =
typename local_matrix_type::index_type::non_const_type;
209 using values_type =
typename local_matrix_type::values_type::non_const_type;
212 RCP<Teuchos::FancyOStream> out;
213 if(
const char* dbg = std::getenv(
"MUELU_REGIONRFACTORY_DEBUG")) {
214 out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
215 out->setShowAllFrontMatter(
false).setShowProcRank(
true);
217 out = Teuchos::getFancyOStream(rcp(
new Teuchos::oblackholestream()));
222 for(
int dim = 0; dim < numDimensions; ++dim) {
223 lCoarseNodesPerDim[dim] = lFineNodesPerDim[dim] / 3 + 1;
225 *out <<
"lCoarseNodesPerDim " << lCoarseNodesPerDim << std::endl;
228 const LO blkSize = A->GetFixedBlockSize();
229 *out <<
"blkSize " << blkSize << std::endl;
233 const LO numRows = blkSize*lCoarseNodesPerDim[0]*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[2];
234 const LO numCols = blkSize*lFineNodesPerDim[0]*lFineNodesPerDim[1]*lFineNodesPerDim[2];
239 RCP<Map> rowMap = MapFactory::Build(A->getRowMap()->lib(),
240 Teuchos::OrdinalTraits<GO>::invalid(),
242 A->getRowMap()->getIndexBase(),
243 A->getRowMap()->getComm());
245 coarseCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(rowMap,
247 Array<ArrayRCP<const real_type> > fineCoordData(numDimensions);
248 Array<ArrayRCP<real_type> > coarseCoordData(numDimensions);
249 for(
int dim = 0; dim < numDimensions; ++dim) {
250 fineCoordData[dim] = fineCoordinates->getData(dim);
251 coarseCoordData[dim] = coarseCoordinates->getDataNonConst(dim);
258 const LO cornerStencilLength = 27;
259 const LO edgeStencilLength = 45;
260 const LO faceStencilLength = 75;
261 const LO interiorStencilLength = 125;
264 const LO numCorners = 8;
265 const LO numEdges = 4*(lCoarseNodesPerDim[0] - 2)
266 + 4*(lCoarseNodesPerDim[1] - 2)
267 + 4*(lCoarseNodesPerDim[2] - 2);
268 const LO numFaces = 2*(lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[1] - 2)
269 + 2*(lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[2] - 2)
270 + 2*(lCoarseNodesPerDim[1] - 2)*(lCoarseNodesPerDim[2] - 2);
271 const LO numInteriors = (lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[1] - 2)
272 *(lCoarseNodesPerDim[2] - 2);
274 const LO nnz = (numCorners*cornerStencilLength + numEdges*edgeStencilLength
275 + numFaces*faceStencilLength + numInteriors*interiorStencilLength)*blkSize;
280 *out <<
"R statistics:" << std::endl
281 <<
" -numRows= " << numRows << std::endl
282 <<
" -numCols= " << numCols << std::endl
283 <<
" -nnz= " << nnz << std::endl;
285 row_map_type row_map(Kokkos::ViewAllocateWithoutInitializing(
"row_map"), numRows + 1);
286 typename row_map_type::HostMirror row_map_h = Kokkos::create_mirror_view(row_map);
288 entries_type entries(Kokkos::ViewAllocateWithoutInitializing(
"entries"), nnz);
289 typename entries_type::HostMirror entries_h = Kokkos::create_mirror_view(entries);
291 values_type values(Kokkos::ViewAllocateWithoutInitializing(
"values"), nnz);
292 typename values_type::HostMirror values_h = Kokkos::create_mirror_view(values);
297 Array<SC> coeffs({1.0/3.0, 2.0/3.0, 1.0, 2.0/3.0, 1.0/3.0});
302 const LO edgeLineOffset = 2*cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength;
303 const LO faceLineOffset = 2*edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength;
304 const LO interiorLineOffset = 2*faceStencilLength
305 + (lCoarseNodesPerDim[0] - 2)*interiorStencilLength;
307 const LO facePlaneOffset = 2*edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset;
308 const LO interiorPlaneOffset = 2*faceLineOffset + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset;
315 LO coordRowIdx = 0, rowIdx = 0, coordColumnOffset = 0, columnOffset = 0, entryOffset = 0;
316 for(LO l = 0; l < blkSize; ++l) {
317 for(LO k = 0; k < 3; ++k) {
318 for(LO j = 0; j < 3; ++j) {
319 for(LO i = 0; i < 3; ++i) {
320 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
321 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
322 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i + 2];
327 for(LO l = 0; l < blkSize; ++l) {
328 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
330 for(
int dim = 0; dim <numDimensions; ++dim) {
331 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
335 coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
336 rowIdx = coordRowIdx*blkSize;
337 coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
338 columnOffset = coordColumnOffset*blkSize;
339 entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
340 for(LO l = 0; l < blkSize; ++l) {
341 for(LO k = 0; k < 3; ++k) {
342 for(LO j = 0; j < 3; ++j) {
343 for(LO i = 0; i < 3; ++i) {
344 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
345 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
346 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i + 2];
351 for(LO l = 0; l < blkSize; ++l) {
352 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
354 for(
int dim = 0; dim <numDimensions; ++dim) {
355 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
359 coordRowIdx = (lCoarseNodesPerDim[0] - 1);
360 rowIdx = coordRowIdx*blkSize;
361 coordColumnOffset = (lFineNodesPerDim[0] - 1);
362 columnOffset = coordColumnOffset*blkSize;
363 entryOffset = (cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength)*blkSize;
364 for(LO l = 0; l < blkSize; ++l) {
365 for(LO k = 0; k < 3; ++k) {
366 for(LO j = 0; j < 3; ++j) {
367 for(LO i = 0; i < 3; ++i) {
368 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
369 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
370 + j*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
371 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i];
376 for(LO l = 0; l < blkSize; ++l) {
377 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
379 for(
int dim = 0; dim <numDimensions; ++dim) {
380 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
384 coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
385 rowIdx = coordRowIdx*blkSize;
386 coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
387 columnOffset = coordColumnOffset*blkSize;
388 entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
389 for(LO l = 0; l < blkSize; ++l) {
390 for(LO k = 0; k < 3; ++k) {
391 for(LO j = 0; j < 3; ++j) {
392 for(LO i = 0; i < 3; ++i) {
393 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
394 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
395 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
400 for(LO l = 0; l < blkSize; ++l) {
401 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
403 for(
int dim = 0; dim <numDimensions; ++dim) {
404 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
408 coordRowIdx = (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0];
409 rowIdx = coordRowIdx*blkSize;
410 coordColumnOffset = (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0];
411 columnOffset = coordColumnOffset*blkSize;
412 entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset)*blkSize;
413 for(LO l = 0; l < blkSize; ++l) {
414 for(LO k = 0; k < 3; ++k) {
415 for(LO j = 0; j < 3; ++j) {
416 for(LO i = 0; i < 3; ++i) {
417 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
418 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
419 + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
420 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i + 2];
425 for(LO l = 0; l < blkSize; ++l) {
426 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
428 for(
int dim = 0; dim <numDimensions; ++dim) {
429 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
433 coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
434 rowIdx = coordRowIdx*blkSize;
435 coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
436 columnOffset = coordColumnOffset*blkSize;
437 entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
438 for(LO l = 0; l < blkSize; ++l) {
439 for(LO k = 0; k < 3; ++k) {
440 for(LO j = 0; j < 3; ++j) {
441 for(LO i = 0; i < 3; ++i) {
442 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
443 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
444 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
449 for(LO l = 0; l < blkSize; ++l) {
450 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
452 for(
int dim = 0; dim <numDimensions; ++dim) {
453 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
457 coordRowIdx = (lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] - 1);
458 rowIdx = coordRowIdx*blkSize;
459 coordColumnOffset = (lFineNodesPerDim[1]*lFineNodesPerDim[0] - 1);
460 columnOffset = coordColumnOffset*blkSize;
461 entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset +
462 cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength)*blkSize;
463 for(LO l = 0; l < blkSize; ++l) {
464 for(LO k = 0; k < 3; ++k) {
465 for(LO j = 0; j < 3; ++j) {
466 for(LO i = 0; i < 3; ++i) {
467 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
468 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
469 + (j - 2)*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
470 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
475 for(LO l = 0; l < blkSize; ++l) {
476 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
478 for(
int dim = 0; dim <numDimensions; ++dim) {
479 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
483 coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
484 rowIdx = coordRowIdx*blkSize;
485 coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
486 columnOffset = coordColumnOffset*blkSize;
487 entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
488 for(LO l = 0; l < blkSize; ++l) {
489 for(LO k = 0; k < 3; ++k) {
490 for(LO j = 0; j < 3; ++j) {
491 for(LO i = 0; i < 3; ++i) {
492 entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
493 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
494 values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
499 for(LO l = 0; l < blkSize; ++l) {
500 row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
502 for(
int dim = 0; dim <numDimensions; ++dim) {
503 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
508 if(lCoarseNodesPerDim[0] - 2 > 0) {
510 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
511 for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[0] - 2; ++edgeIdx) {
514 coordRowIdx = (edgeIdx + 1);
515 rowIdx = coordRowIdx*blkSize;
516 coordColumnOffset = (edgeIdx + 1)*3;
517 columnOffset = coordColumnOffset*blkSize;
518 entryOffset = (cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
519 for(LO l = 0; l < blkSize; ++l) {
520 for(LO k = 0; k < 3; ++k) {
521 for(LO j = 0; j < 3; ++j) {
522 for(LO i = 0; i < 5; ++i) {
523 entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
524 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
525 values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i];
530 for(LO l = 0; l < blkSize; ++l) {
531 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
533 for(
int dim = 0; dim <numDimensions; ++dim) {
534 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
538 coordRowIdx = ((lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0] + edgeIdx + 1);
539 rowIdx = coordRowIdx*blkSize;
540 coordColumnOffset = ((lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0] + (edgeIdx + 1)*3);
541 columnOffset = coordColumnOffset*blkSize;
542 entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset
543 + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
544 for(LO l = 0; l < blkSize; ++l) {
545 for(LO k = 0; k < 3; ++k) {
546 for(LO j = 0; j < 3; ++j) {
547 for(LO i = 0; i < 5; ++i) {
548 entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
549 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
550 values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
555 for(LO l = 0; l < blkSize; ++l) {
556 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
558 for(
int dim = 0; dim <numDimensions; ++dim) {
559 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
563 coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
565 rowIdx = coordRowIdx*blkSize;
566 coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
568 columnOffset = coordColumnOffset*blkSize;
569 entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
570 + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
571 for(LO l = 0; l < blkSize; ++l) {
572 for(LO k = 0; k < 3; ++k) {
573 for(LO j = 0; j < 3; ++j) {
574 for(LO i = 0; i < 5; ++i) {
575 entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
576 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
577 values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
582 for(LO l = 0; l < blkSize; ++l) {
583 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
585 for(
int dim = 0; dim <numDimensions; ++dim) {
586 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
590 coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
591 + (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0] + edgeIdx + 1);
592 rowIdx = coordRowIdx*blkSize;
593 coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
594 + (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0] + (edgeIdx + 1)*3);
595 columnOffset = coordColumnOffset*blkSize;
596 entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
597 + edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset
598 + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
599 for(LO l = 0; l < blkSize; ++l) {
600 for(LO k = 0; k < 3; ++k) {
601 for(LO j = 0; j < 3; ++j) {
602 for(LO i = 0; i < 5; ++i) {
603 entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
604 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
605 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
606 values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
611 for(LO l = 0; l < blkSize; ++l) {
612 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
614 for(
int dim = 0; dim <numDimensions; ++dim) {
615 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
621 if(lCoarseNodesPerDim[1] - 2 > 0) {
623 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
624 for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[1] - 2; ++edgeIdx) {
627 coordRowIdx = (edgeIdx + 1)*lCoarseNodesPerDim[0];
628 rowIdx = coordRowIdx*blkSize;
629 coordColumnOffset = (edgeIdx + 1)*3*lFineNodesPerDim[0];
630 columnOffset = coordColumnOffset*blkSize;
631 entryOffset = (edgeLineOffset + edgeIdx*faceLineOffset)*blkSize;
632 for(LO l = 0; l < blkSize; ++l) {
633 for(LO k = 0; k < 3; ++k) {
634 for(LO j = 0; j < 5; ++j) {
635 for(LO i = 0; i < 3; ++i) {
636 entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
637 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
638 values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i + 2];
643 for(LO l = 0; l < blkSize; ++l) {
644 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
646 for(
int dim = 0; dim <numDimensions; ++dim) {
647 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
651 coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
652 rowIdx = coordRowIdx*blkSize;
653 coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
654 columnOffset = coordColumnOffset*blkSize;
655 entryOffset = (edgeLineOffset + edgeIdx*faceLineOffset
656 + edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength)*blkSize;
657 for(LO l = 0; l < blkSize; ++l) {
658 for(LO k = 0; k < 3; ++k) {
659 for(LO j = 0; j < 5; ++j) {
660 for(LO i = 0; i < 3; ++i) {
661 entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
662 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
663 values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
668 for(LO l = 0; l < blkSize; ++l) {
669 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
671 for(
int dim = 0; dim <numDimensions; ++dim) {
672 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
676 coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
677 + (edgeIdx + 1)*lCoarseNodesPerDim[0]);
678 rowIdx = coordRowIdx*blkSize;
679 coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
680 + (edgeIdx + 1)*3*lFineNodesPerDim[0]);
681 columnOffset = coordColumnOffset*blkSize;
682 entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
683 + edgeLineOffset + edgeIdx*faceLineOffset)*blkSize;
684 for(LO l = 0; l < blkSize; ++l) {
685 for(LO k = 0; k < 3; ++k) {
686 for(LO j = 0; j < 5; ++j) {
687 for(LO i = 0; i < 3; ++i) {
688 entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
689 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
690 values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
695 for(LO l = 0; l < blkSize; ++l) {
696 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
698 for(
int dim = 0; dim <numDimensions; ++dim) {
699 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
703 coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
704 + (edgeIdx + 1)*lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
705 rowIdx = coordRowIdx*blkSize;
706 coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
707 + (edgeIdx + 1)*3*lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
708 columnOffset = coordColumnOffset*blkSize;
709 entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
710 + edgeLineOffset + edgeIdx*faceLineOffset
711 + edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength)*blkSize;
712 for(LO l = 0; l < blkSize; ++l) {
713 for(LO k = 0; k < 3; ++k) {
714 for(LO j = 0; j < 5; ++j) {
715 for(LO i = 0; i < 3; ++i) {
716 entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
717 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
718 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
719 values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
724 for(LO l = 0; l < blkSize; ++l) {
725 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
727 for(
int dim = 0; dim <numDimensions; ++dim) {
728 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
734 if(lCoarseNodesPerDim[2] - 2 > 0) {
736 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
737 for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[2] - 2; ++edgeIdx) {
740 coordRowIdx = (edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
741 rowIdx = coordRowIdx*blkSize;
742 coordColumnOffset = (edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0];
743 columnOffset = coordColumnOffset*blkSize;
744 entryOffset = (facePlaneOffset + edgeIdx*interiorPlaneOffset)*blkSize;
745 for(LO l = 0; l < blkSize; ++l) {
746 for(LO k = 0; k < 5; ++k) {
747 for(LO j = 0; j < 3; ++j) {
748 for(LO i = 0; i < 3; ++i) {
749 entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
750 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
751 values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i + 2];
756 for(LO l = 0; l < blkSize; ++l) {
757 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
759 for(
int dim = 0; dim <numDimensions; ++dim) {
760 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
764 coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
765 + lCoarseNodesPerDim[0] - 1);
766 rowIdx = coordRowIdx*blkSize;
767 coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
768 + lFineNodesPerDim[0] - 1);
769 columnOffset = coordColumnOffset*blkSize;
770 entryOffset = (facePlaneOffset + faceLineOffset - edgeStencilLength
771 + edgeIdx*interiorPlaneOffset)*blkSize;
772 for(LO l = 0; l < blkSize; ++l) {
773 for(LO k = 0; k < 5; ++k) {
774 for(LO j = 0; j < 3; ++j) {
775 for(LO i = 0; i < 3; ++i) {
776 entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
777 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
778 values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
783 for(LO l = 0; l < blkSize; ++l) {
784 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
786 for(
int dim = 0; dim <numDimensions; ++dim) {
787 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
791 coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
792 + (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0]);
793 rowIdx = coordRowIdx*blkSize;
794 coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
795 + (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0]);
796 columnOffset = coordColumnOffset*blkSize;
797 entryOffset = (facePlaneOffset + edgeIdx*interiorPlaneOffset + faceLineOffset
798 + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset)*blkSize;
799 for(LO l = 0; l < blkSize; ++l) {
800 for(LO k = 0; k < 5; ++k) {
801 for(LO j = 0; j < 3; ++j) {
802 for(LO i = 0; i < 3; ++i) {
803 entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
804 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
805 values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
810 for(LO l = 0; l < blkSize; ++l) {
811 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
813 for(
int dim = 0; dim <numDimensions; ++dim) {
814 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
818 coordRowIdx = ((edgeIdx + 2)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] - 1);
819 rowIdx = coordRowIdx*blkSize;
820 coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
821 + lFineNodesPerDim[1]*lFineNodesPerDim[0] - 1);
822 columnOffset = coordColumnOffset*blkSize;
823 entryOffset = (facePlaneOffset + (edgeIdx + 1)*interiorPlaneOffset - edgeStencilLength)*blkSize;
824 for(LO l = 0; l < blkSize; ++l) {
825 for(LO k = 0; k < 5; ++k) {
826 for(LO j = 0; j < 3; ++j) {
827 for(LO i = 0; i < 3; ++i) {
828 entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
829 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
830 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
831 values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
836 for(LO l = 0; l < blkSize; ++l) {
837 row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
839 for(
int dim = 0; dim <numDimensions; ++dim) {
840 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
846 if((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[1] - 2 > 0)) {
848 Array<LO> gridIdx(3);
849 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
850 for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[1]-2)*(lCoarseNodesPerDim[0]-2); ++faceIdx) {
853 coordRowIdx = ((gridIdx[1] + 1)*lCoarseNodesPerDim[0] + gridIdx[0] + 1);
854 rowIdx = coordRowIdx*blkSize;
855 coordColumnOffset = 3*((gridIdx[1] + 1)*lFineNodesPerDim[0] + gridIdx[0] + 1);
856 columnOffset = coordColumnOffset*blkSize;
857 entryOffset = (edgeLineOffset + edgeStencilLength
858 + gridIdx[1]*faceLineOffset + gridIdx[0]*faceStencilLength)*blkSize;
859 for(LO l = 0; l < blkSize; ++l) {
860 for(LO k = 0; k < 3; ++k) {
861 for(LO j = 0; j < 5; ++j) {
862 for(LO i = 0; i < 5; ++i) {
863 entries_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = columnOffset
864 + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
865 values_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
870 for(LO l = 0; l < blkSize; ++l) {
871 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
873 for(
int dim = 0; dim <numDimensions; ++dim) {
874 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
878 coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
879 rowIdx = coordRowIdx*blkSize;
880 coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
881 columnOffset = coordColumnOffset*blkSize;
882 entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
883 for(LO l = 0; l < blkSize; ++l) {
884 for(LO k = 0; k < 3; ++k) {
885 for(LO j = 0; j < 5; ++j) {
886 for(LO i = 0; i < 5; ++i) {
887 entries_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = columnOffset
888 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
889 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
890 values_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
895 for(LO l = 0; l < blkSize; ++l) {
896 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
898 for(
int dim = 0; dim <numDimensions; ++dim) {
899 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
906 if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
914 if((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
916 Array<LO> gridIdx(3);
917 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
918 for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[2]-2)*(lCoarseNodesPerDim[0]-2); ++faceIdx) {
921 coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] + (gridIdx[0] + 1));
922 rowIdx = coordRowIdx*blkSize;
923 coordColumnOffset = ((gridIdx[2] + 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
925 columnOffset = coordColumnOffset*blkSize;
926 entryOffset = (facePlaneOffset + gridIdx[2]*interiorPlaneOffset + edgeStencilLength
927 + gridIdx[0]*faceStencilLength)*blkSize;
928 for(LO l = 0; l < blkSize; ++l) {
929 for(LO k = 0; k < 5; ++k) {
930 for(LO j = 0; j < 3; ++j) {
931 for(LO i = 0; i < 5; ++i) {
932 entries_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = columnOffset
933 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
934 values_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
939 for(LO l = 0; l < blkSize; ++l) {
940 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
942 for(
int dim = 0; dim <numDimensions; ++dim) {
943 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
947 coordRowIdx += (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0];
948 rowIdx = coordRowIdx*blkSize;
949 coordColumnOffset += (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0];
950 columnOffset = coordColumnOffset*blkSize;
951 entryOffset += (faceLineOffset + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset)*blkSize;
952 for(LO l = 0; l < blkSize; ++l) {
953 for(LO k = 0; k < 5; ++k) {
954 for(LO j = 0; j < 3; ++j) {
955 for(LO i = 0; i < 5; ++i) {
956 entries_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = columnOffset
957 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
958 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
959 values_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
964 for(LO l = 0; l < blkSize; ++l) {
965 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
967 for(
int dim = 0; dim <numDimensions; ++dim) {
968 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
975 if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
983 if((lCoarseNodesPerDim[1] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
985 Array<LO> gridIdx(3);
986 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
987 for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[2]-2)*(lCoarseNodesPerDim[1]-2); ++faceIdx) {
990 coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
991 + (gridIdx[1] + 1)*lCoarseNodesPerDim[0]);
992 rowIdx = coordRowIdx*blkSize;
993 coordColumnOffset = ((gridIdx[2] + 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
994 + (gridIdx[1] + 1)*lFineNodesPerDim[0])*3;
995 columnOffset = coordColumnOffset*blkSize;
996 entryOffset = (facePlaneOffset + gridIdx[2]*interiorPlaneOffset + faceLineOffset
997 + gridIdx[1]*interiorLineOffset)*blkSize;
998 for(LO l = 0; l < blkSize; ++l) {
999 for(LO k = 0; k < 5; ++k) {
1000 for(LO j = 0; j < 5; ++j) {
1001 for(LO i = 0; i < 3; ++i) {
1002 entries_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = columnOffset
1003 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
1004 values_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
1009 for(LO l = 0; l < blkSize; ++l) {
1010 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
1012 for(
int dim = 0; dim <numDimensions; ++dim) {
1013 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1017 coordRowIdx += (lCoarseNodesPerDim[0] - 1);
1018 rowIdx = coordRowIdx*blkSize;
1019 coordColumnOffset += (lFineNodesPerDim[0] - 1);
1020 columnOffset = coordColumnOffset*blkSize;
1021 entryOffset += (faceStencilLength + (lCoarseNodesPerDim[0] - 2)*interiorStencilLength)*blkSize;
1022 for(LO l = 0; l < blkSize; ++l) {
1023 for(LO k = 0; k < 5; ++k) {
1024 for(LO j = 0; j < 5; ++j) {
1025 for(LO i = 0; i < 3; ++i) {
1026 entries_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = columnOffset
1027 + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1028 + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
1029 values_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
1034 for(LO l = 0; l < blkSize; ++l) {
1035 row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
1037 for(
int dim = 0; dim <numDimensions; ++dim) {
1038 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1045 if(gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
1052 if(numInteriors > 0) {
1057 LO countRowEntries = 0;
1058 Array<LO> coordColumnOffsets(125);
1059 for(LO k = -2; k < 3; ++k) {
1060 for(LO j = -2; j < 3; ++j) {
1061 for(LO i = -2; i < 3; ++i) {
1062 coordColumnOffsets[countRowEntries] = k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1063 + j*lFineNodesPerDim[0] + i;
1070 Array<SC> interiorValues(125);
1071 for(LO k = 0; k < 5; ++k) {
1072 for(LO j = 0; j < 5; ++j) {
1073 for(LO i = 0; i < 5; ++i) {
1074 interiorValues[countValues] = coeffs[k]*coeffs[j]*coeffs[i];
1080 LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
1081 Array<LO> gridIdx(3);
1082 for(LO interiorIdx = 0; interiorIdx < numInteriors; ++interiorIdx) {
1083 coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[0]*lCoarseNodesPerDim[1]
1084 + (gridIdx[1] + 1)*lCoarseNodesPerDim[0]
1086 rowIdx = coordRowIdx*blkSize;
1087 coordColumnOffset = ((gridIdx[2] + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1088 + (gridIdx[1] + 1)*3*lFineNodesPerDim[0] + (gridIdx[0] + 1)*3);
1089 columnOffset = coordColumnOffset*blkSize;
1090 entryOffset = (facePlaneOffset + faceLineOffset + faceStencilLength
1091 + gridIdx[2]*interiorPlaneOffset + gridIdx[1]*interiorLineOffset
1092 + gridIdx[0]*interiorStencilLength)*blkSize;
1093 for(LO l = 0; l < blkSize; ++l) {
1094 row_map_h(rowIdx + 1 + l) = entryOffset + interiorStencilLength*(l+1);
1099 for(LO l = 0; l < blkSize; ++l) {
1100 for(LO entryIdx = 0; entryIdx < interiorStencilLength; ++entryIdx) {
1101 entries_h(entryOffset + entryIdx + interiorStencilLength*l) = columnOffset + coordColumnOffsets[entryIdx]*blkSize + l;
1102 values_h(entryOffset + entryIdx + interiorStencilLength*l) = interiorValues[entryIdx];
1105 for(
int dim = 0; dim <numDimensions; ++dim) {
1106 coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1113 if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
1116 if(gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
1124 Kokkos::deep_copy(row_map, row_map_h);
1125 Kokkos::deep_copy(entries, entries_h);
1126 Kokkos::deep_copy(values, values_h);
1128 local_graph_type localGraph(entries, row_map);
1129 local_matrix_type localR(
"R", numCols, values, localGraph);
1131 R = MatrixFactory::Build(localR,
1142 #define MUELU_REGIONRFACTORY_SHORT 1143 #endif //ifdef HAVE_MUELU_KOKKOS_REFACTOR 1144 #endif // MUELU_REGIONRFACTORY_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Namespace for MueLu classes and methods.
static const NoFactory * get()
MueLu::DefaultGlobalOrdinal GlobalOrdinal
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)