MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_RegionRFactory_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_REGIONRFACTORY_DEF_HPP
11 #define MUELU_REGIONRFACTORY_DEF_HPP
12 
13 #include <Xpetra_Matrix.hpp>
14 #include <Xpetra_CrsGraphFactory.hpp>
15 
16 // #include "MueLu_PFactory.hpp"
17 // #include "MueLu_FactoryManagerBase.hpp"
18 
19 #include "MueLu_RegionRFactory_decl.hpp"
20 
21 namespace MueLu {
22 
23 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
24 RCP<const ParameterList> RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
25  GetValidParameterList() const {
26  RCP<ParameterList> validParamList = rcp(new ParameterList());
27 
28  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
29  "Generating factory of the matrix A");
30  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
31  "Number of spatial dimensions in the problem.");
32  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
33  "Number of local nodes per spatial dimension on the fine grid.");
34  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
35  "Fine level nullspace used to construct the coarse level nullspace.");
36  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
37  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
38  validParamList->set<bool>("keep coarse coords", false, "Flag to keep coordinates for special coarse grid solve");
39 
40  return validParamList;
41 } // GetValidParameterList()
42 
43 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
44 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
45  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
46  Input(fineLevel, "A");
47  Input(fineLevel, "numDimensions");
48  Input(fineLevel, "lNodesPerDim");
49  Input(fineLevel, "Nullspace");
50  Input(fineLevel, "Coordinates");
51 
52 } // DeclareInput()
53 
54 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
56  Build(Level& fineLevel, Level& coarseLevel) const {
57  // Set debug outputs based on environment variable
58  RCP<Teuchos::FancyOStream> out;
59  if (const char* dbg = std::getenv("MUELU_REGIONRFACTORY_DEBUG")) {
60  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
61  out->setShowAllFrontMatter(false).setShowProcRank(true);
62  } else {
63  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
64  }
65 
66  *out << "Starting RegionRFactory::Build." << std::endl;
67 
68  // First get the inputs from the fineLevel
69  const int numDimensions = Get<int>(fineLevel, "numDimensions");
70  Array<LO> lFineNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
71  {
72  Array<LO> lNodesPerDim = Get<Array<LO> >(fineLevel, "lNodesPerDim");
73  for (int dim = 0; dim < numDimensions; ++dim) {
74  lFineNodesPerDim[dim] = lNodesPerDim[dim];
75  }
76  }
77  *out << "numDimensions " << numDimensions << " and lFineNodesPerDim: " << lFineNodesPerDim
78  << std::endl;
79 
80  // Let us check that the inputs verify our assumptions
81  if (numDimensions < 1 || numDimensions > 3) {
82  throw std::runtime_error("numDimensions must be 1, 2 or 3!");
83  }
84  for (int dim = 0; dim < numDimensions; ++dim) {
85  if (lFineNodesPerDim[dim] % 3 != 1) {
86  throw std::runtime_error("The number of fine node in each direction need to be 3n+1");
87  }
88  }
89  Array<LO> lCoarseNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
90 
91  const RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
92 
93  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
94  fineCoordinates = Get<RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
95  if (static_cast<int>(fineCoordinates->getNumVectors()) != numDimensions) {
96  throw std::runtime_error("The number of vectors in the coordinates is not equal to numDimensions!");
97  }
98 
99  // Let us create R and pass it down to the
100  // appropriate specialization and see what we
101  // get back!
102  RCP<Matrix> R;
103 
104  if (numDimensions == 1) {
105  throw std::runtime_error("RegionRFactory no implemented for 1D case yet.");
106  } else if (numDimensions == 2) {
107  throw std::runtime_error("RegionRFactory no implemented for 2D case yet.");
108  } else if (numDimensions == 3) {
109  Build3D(numDimensions, lFineNodesPerDim, A, fineCoordinates,
110  R, coarseCoordinates, lCoarseNodesPerDim);
111  }
112 
113  const Teuchos::ParameterList& pL = GetParameterList();
114 
115  // Reuse pattern if available (multiple solve)
116  RCP<ParameterList> Tparams;
117  if (pL.isSublist("matrixmatrix: kernel params"))
118  Tparams = rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
119  else
120  Tparams = rcp(new ParameterList);
121 
122  // R->describe(*out, Teuchos::VERB_EXTREME);
123  *out << "Compute P=R^t" << std::endl;
124  // By default, we don't need global constants for transpose
125  Tparams->set("compute global constants: temporaries", Tparams->get("compute global constants: temporaries", false));
126  Tparams->set("compute global constants", Tparams->get("compute global constants", false));
127  std::string label = "MueLu::RegionR-transR" + Teuchos::toString(coarseLevel.GetLevelID());
128  RCP<Matrix> P = Utilities::Transpose(*R, true, label, Tparams);
129 
130  *out << "Compute coarse nullspace" << std::endl;
131  RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel, "Nullspace");
132  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(R->getRowMap(),
133  fineNullspace->getNumVectors());
134  R->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
135  Teuchos::ScalarTraits<SC>::zero());
136 
137  *out << "Set data on coarse level" << std::endl;
138  Set(coarseLevel, "numDimensions", numDimensions);
139  Set(coarseLevel, "lNodesPerDim", lCoarseNodesPerDim);
140  Set(coarseLevel, "Nullspace", coarseNullspace);
141  Set(coarseLevel, "Coordinates", coarseCoordinates);
142  if (pL.get<bool>("keep coarse coords")) {
143  coarseLevel.Set<RCP<realvaluedmultivector_type> >("Coordinates2", coarseCoordinates, NoFactory::get());
144  }
145 
146  R->SetFixedBlockSize(A->GetFixedBlockSize());
147  P->SetFixedBlockSize(A->GetFixedBlockSize());
148 
149  Set(coarseLevel, "R", R);
150  Set(coarseLevel, "P", P);
151 
152 } // Build()
153 
154 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
155 void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
156  Build3D(const int numDimensions,
157  Teuchos::Array<LocalOrdinal>& lFineNodesPerDim,
158  const RCP<Matrix>& A,
159  const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> >& fineCoordinates,
160  RCP<Matrix>& R,
161  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> >& coarseCoordinates,
162  Teuchos::Array<LocalOrdinal>& lCoarseNodesPerDim) const {
163  using local_matrix_type = typename CrsMatrix::local_matrix_type;
164  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
165  using row_map_type = typename local_matrix_type::row_map_type::non_const_type;
166  using entries_type = typename local_matrix_type::index_type::non_const_type;
167  using values_type = typename local_matrix_type::values_type::non_const_type;
168 
169  // Set debug outputs based on environment variable
170  RCP<Teuchos::FancyOStream> out;
171  if (const char* dbg = std::getenv("MUELU_REGIONRFACTORY_DEBUG")) {
172  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
173  out->setShowAllFrontMatter(false).setShowProcRank(true);
174  } else {
175  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
176  }
177 
178  // Now compute number of coarse grid points
179  for (int dim = 0; dim < numDimensions; ++dim) {
180  lCoarseNodesPerDim[dim] = lFineNodesPerDim[dim] / 3 + 1;
181  }
182  *out << "lCoarseNodesPerDim " << lCoarseNodesPerDim << std::endl;
183 
184  // Grab the block size here and multiply all existing offsets by it
185  const LO blkSize = A->GetFixedBlockSize();
186  *out << "blkSize " << blkSize << std::endl;
187 
188  // Based on lCoarseNodesPerDim and lFineNodesPerDim
189  // we can compute numRows, numCols and NNZ for R
190  const LO numRows = blkSize * lCoarseNodesPerDim[0] * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[2];
191  const LO numCols = blkSize * lFineNodesPerDim[0] * lFineNodesPerDim[1] * lFineNodesPerDim[2];
192 
193  // Create the coarse coordinates multivector
194  // so we can fill it on the fly while computing
195  // the restriction operator
196  RCP<Map> rowMap = MapFactory::Build(A->getRowMap()->lib(),
197  Teuchos::OrdinalTraits<GO>::invalid(),
198  numRows,
199  A->getRowMap()->getIndexBase(),
200  A->getRowMap()->getComm());
201 
202  RCP<Map> coordRowMap = MapFactory::Build(A->getRowMap()->lib(),
203  Teuchos::OrdinalTraits<GO>::invalid(),
204  lCoarseNodesPerDim[0] * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[2],
205  A->getRowMap()->getIndexBase(),
206  A->getRowMap()->getComm());
207 
208  coarseCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(coordRowMap,
209  numDimensions);
210  Array<ArrayRCP<const real_type> > fineCoordData(numDimensions);
211  Array<ArrayRCP<real_type> > coarseCoordData(numDimensions);
212  for (int dim = 0; dim < numDimensions; ++dim) {
213  fineCoordData[dim] = fineCoordinates->getData(dim);
214  coarseCoordData[dim] = coarseCoordinates->getDataNonConst(dim);
215  }
216 
217  // Let us set some parameter that will be useful
218  // while constructing R
219 
220  // Length of interpolation stencils based on geometry
221  const LO cornerStencilLength = 27;
222  const LO edgeStencilLength = 45;
223  const LO faceStencilLength = 75;
224  const LO interiorStencilLength = 125;
225 
226  // Number of corner, edge, face and interior nodes
227  const LO numCorners = 8;
228  const LO numEdges = 4 * (lCoarseNodesPerDim[0] - 2) + 4 * (lCoarseNodesPerDim[1] - 2) + 4 * (lCoarseNodesPerDim[2] - 2);
229  const LO numFaces = 2 * (lCoarseNodesPerDim[0] - 2) * (lCoarseNodesPerDim[1] - 2) + 2 * (lCoarseNodesPerDim[0] - 2) * (lCoarseNodesPerDim[2] - 2) + 2 * (lCoarseNodesPerDim[1] - 2) * (lCoarseNodesPerDim[2] - 2);
230  const LO numInteriors = (lCoarseNodesPerDim[0] - 2) * (lCoarseNodesPerDim[1] - 2) * (lCoarseNodesPerDim[2] - 2);
231 
232  const LO nnz = (numCorners * cornerStencilLength + numEdges * edgeStencilLength + numFaces * faceStencilLength + numInteriors * interiorStencilLength) * blkSize;
233 
234  // Having the number of rows and columns we can genrate
235  // the appropriate maps for our operator.
236 
237  *out << "R statistics:" << std::endl
238  << " -numRows= " << numRows << std::endl
239  << " -numCols= " << numCols << std::endl
240  << " -nnz= " << nnz << std::endl;
241 
242  row_map_type row_map(Kokkos::ViewAllocateWithoutInitializing("row_map"), numRows + 1);
243  typename row_map_type::HostMirror row_map_h = Kokkos::create_mirror_view(row_map);
244 
245  entries_type entries(Kokkos::ViewAllocateWithoutInitializing("entries"), nnz);
246  typename entries_type::HostMirror entries_h = Kokkos::create_mirror_view(entries);
247 
248  values_type values(Kokkos::ViewAllocateWithoutInitializing("values"), nnz);
249  typename values_type::HostMirror values_h = Kokkos::create_mirror_view(values);
250 
251  // Compute the basic interpolation
252  // coefficients for 1D rate of 3
253  // coarsening.
254  Array<SC> coeffs({1.0 / 3.0, 2.0 / 3.0, 1.0, 2.0 / 3.0, 1.0 / 3.0});
255  row_map_h(0) = 0;
256 
257  // Define some offsets that
258  // will be needed often later on
259  const LO edgeLineOffset = 2 * cornerStencilLength + (lCoarseNodesPerDim[0] - 2) * edgeStencilLength;
260  const LO faceLineOffset = 2 * edgeStencilLength + (lCoarseNodesPerDim[0] - 2) * faceStencilLength;
261  const LO interiorLineOffset = 2 * faceStencilLength + (lCoarseNodesPerDim[0] - 2) * interiorStencilLength;
262 
263  const LO facePlaneOffset = 2 * edgeLineOffset + (lCoarseNodesPerDim[1] - 2) * faceLineOffset;
264  const LO interiorPlaneOffset = 2 * faceLineOffset + (lCoarseNodesPerDim[1] - 2) * interiorLineOffset;
265 
266  // Let us take care of the corners
267  // first since we always have
268  // corners to deal with!
269  {
270  // Corner 1
271  LO coordRowIdx = 0, rowIdx = 0, coordColumnOffset = 0, columnOffset = 0, entryOffset = 0;
272  for (LO l = 0; l < blkSize; ++l) {
273  for (LO k = 0; k < 3; ++k) {
274  for (LO j = 0; j < 3; ++j) {
275  for (LO i = 0; i < 3; ++i) {
276  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i) * blkSize + l;
277  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k + 2] * coeffs[j + 2] * coeffs[i + 2];
278  }
279  }
280  }
281  }
282  for (LO l = 0; l < blkSize; ++l) {
283  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
284  }
285  for (int dim = 0; dim < numDimensions; ++dim) {
286  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
287  }
288 
289  // Corner 5
290  coordRowIdx += (lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
291  rowIdx = coordRowIdx * blkSize;
292  coordColumnOffset += (lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0];
293  columnOffset = coordColumnOffset * blkSize;
294  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset) * blkSize;
295  for (LO l = 0; l < blkSize; ++l) {
296  for (LO k = 0; k < 3; ++k) {
297  for (LO j = 0; j < 3; ++j) {
298  for (LO i = 0; i < 3; ++i) {
299  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i) * blkSize + l;
300  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i + 2];
301  }
302  }
303  }
304  }
305  for (LO l = 0; l < blkSize; ++l) {
306  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
307  }
308  for (int dim = 0; dim < numDimensions; ++dim) {
309  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
310  }
311 
312  // Corner 2
313  coordRowIdx = (lCoarseNodesPerDim[0] - 1);
314  rowIdx = coordRowIdx * blkSize;
315  coordColumnOffset = (lFineNodesPerDim[0] - 1);
316  columnOffset = coordColumnOffset * blkSize;
317  entryOffset = (cornerStencilLength + (lCoarseNodesPerDim[0] - 2) * edgeStencilLength) * blkSize;
318  for (LO l = 0; l < blkSize; ++l) {
319  for (LO k = 0; k < 3; ++k) {
320  for (LO j = 0; j < 3; ++j) {
321  for (LO i = 0; i < 3; ++i) {
322  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + (i - 2)) * blkSize + l;
323  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k + 2] * coeffs[j + 2] * coeffs[i];
324  }
325  }
326  }
327  }
328  for (LO l = 0; l < blkSize; ++l) {
329  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
330  }
331  for (int dim = 0; dim < numDimensions; ++dim) {
332  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
333  }
334 
335  // Corner 6
336  coordRowIdx += (lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
337  rowIdx = coordRowIdx * blkSize;
338  coordColumnOffset += (lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0];
339  columnOffset = coordColumnOffset * blkSize;
340  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset) * blkSize;
341  for (LO l = 0; l < blkSize; ++l) {
342  for (LO k = 0; k < 3; ++k) {
343  for (LO j = 0; j < 3; ++j) {
344  for (LO i = 0; i < 3; ++i) {
345  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i - 2) * blkSize + l;
346  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i];
347  }
348  }
349  }
350  }
351  for (LO l = 0; l < blkSize; ++l) {
352  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
353  }
354  for (int dim = 0; dim < numDimensions; ++dim) {
355  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
356  }
357 
358  // Corner 3
359  coordRowIdx = (lCoarseNodesPerDim[1] - 1) * lCoarseNodesPerDim[0];
360  rowIdx = coordRowIdx * blkSize;
361  coordColumnOffset = (lFineNodesPerDim[1] - 1) * lFineNodesPerDim[0];
362  columnOffset = coordColumnOffset * blkSize;
363  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2) * faceLineOffset) * 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 + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
369  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i + 2];
370  }
371  }
372  }
373  }
374  for (LO l = 0; l < blkSize; ++l) {
375  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
376  }
377  for (int dim = 0; dim < numDimensions; ++dim) {
378  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
379  }
380 
381  // Corner 7
382  coordRowIdx += (lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
383  rowIdx = coordRowIdx * blkSize;
384  coordColumnOffset += (lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0];
385  columnOffset = coordColumnOffset * blkSize;
386  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset) * blkSize;
387  for (LO l = 0; l < blkSize; ++l) {
388  for (LO k = 0; k < 3; ++k) {
389  for (LO j = 0; j < 3; ++j) {
390  for (LO i = 0; i < 3; ++i) {
391  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
392  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i + 2];
393  }
394  }
395  }
396  }
397  for (LO l = 0; l < blkSize; ++l) {
398  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
399  }
400  for (int dim = 0; dim < numDimensions; ++dim) {
401  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
402  }
403 
404  // Corner 4
405  coordRowIdx = (lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] - 1);
406  rowIdx = coordRowIdx * blkSize;
407  coordColumnOffset = (lFineNodesPerDim[1] * lFineNodesPerDim[0] - 1);
408  columnOffset = coordColumnOffset * blkSize;
409  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2) * faceLineOffset +
410  cornerStencilLength + (lCoarseNodesPerDim[0] - 2) * edgeStencilLength) *
411  blkSize;
412  for (LO l = 0; l < blkSize; ++l) {
413  for (LO k = 0; k < 3; ++k) {
414  for (LO j = 0; j < 3; ++j) {
415  for (LO i = 0; i < 3; ++i) {
416  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + (i - 2)) * blkSize + l;
417  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i];
418  }
419  }
420  }
421  }
422  for (LO l = 0; l < blkSize; ++l) {
423  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
424  }
425  for (int dim = 0; dim < numDimensions; ++dim) {
426  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
427  }
428 
429  // Corner 8
430  coordRowIdx += (lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
431  rowIdx = coordRowIdx * blkSize;
432  coordColumnOffset += (lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0];
433  columnOffset = coordColumnOffset * blkSize;
434  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset) * blkSize;
435  for (LO l = 0; l < blkSize; ++l) {
436  for (LO k = 0; k < 3; ++k) {
437  for (LO j = 0; j < 3; ++j) {
438  for (LO i = 0; i < 3; ++i) {
439  entries_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + (i - 2)) * blkSize + l;
440  values_h(entryOffset + k * 9 + j * 3 + i + cornerStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
441  }
442  }
443  }
444  }
445  for (LO l = 0; l < blkSize; ++l) {
446  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength * (l + 1);
447  }
448  for (int dim = 0; dim < numDimensions; ++dim) {
449  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
450  }
451  } // Corners are done!
452 
453  // Edges along 0 direction
454  if (lCoarseNodesPerDim[0] - 2 > 0) {
455  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
456  for (LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[0] - 2; ++edgeIdx) {
457  // Edge 0
458  coordRowIdx = (edgeIdx + 1);
459  rowIdx = coordRowIdx * blkSize;
460  coordColumnOffset = (edgeIdx + 1) * 3;
461  columnOffset = coordColumnOffset * blkSize;
462  entryOffset = (cornerStencilLength + edgeIdx * 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 < 5; ++i) {
467  entries_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i - 2) * blkSize + l;
468  values_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = coeffs[k + 2] * coeffs[j + 2] * coeffs[i];
469  }
470  }
471  }
472  }
473  for (LO l = 0; l < blkSize; ++l) {
474  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
475  }
476  for (int dim = 0; dim < numDimensions; ++dim) {
477  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
478  }
479 
480  // Edge 1
481  coordRowIdx = ((lCoarseNodesPerDim[1] - 1) * lCoarseNodesPerDim[0] + edgeIdx + 1);
482  rowIdx = coordRowIdx * blkSize;
483  coordColumnOffset = ((lFineNodesPerDim[1] - 1) * lFineNodesPerDim[0] + (edgeIdx + 1) * 3);
484  columnOffset = coordColumnOffset * blkSize;
485  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2) * faceLineOffset + cornerStencilLength + edgeIdx * edgeStencilLength) * blkSize;
486  for (LO l = 0; l < blkSize; ++l) {
487  for (LO k = 0; k < 3; ++k) {
488  for (LO j = 0; j < 3; ++j) {
489  for (LO i = 0; i < 5; ++i) {
490  entries_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
491  values_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i];
492  }
493  }
494  }
495  }
496  for (LO l = 0; l < blkSize; ++l) {
497  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
498  }
499  for (int dim = 0; dim < numDimensions; ++dim) {
500  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
501  }
502 
503  // Edge 2
504  coordRowIdx = ((lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + edgeIdx + 1);
505  rowIdx = coordRowIdx * blkSize;
506  coordColumnOffset = ((lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (edgeIdx + 1) * 3);
507  columnOffset = coordColumnOffset * blkSize;
508  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset + cornerStencilLength + edgeIdx * edgeStencilLength) * blkSize;
509  for (LO l = 0; l < blkSize; ++l) {
510  for (LO k = 0; k < 3; ++k) {
511  for (LO j = 0; j < 3; ++j) {
512  for (LO i = 0; i < 5; ++i) {
513  entries_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i - 2) * blkSize + l;
514  values_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i];
515  }
516  }
517  }
518  }
519  for (LO l = 0; l < blkSize; ++l) {
520  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
521  }
522  for (int dim = 0; dim < numDimensions; ++dim) {
523  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
524  }
525 
526  // Edge 3
527  coordRowIdx = ((lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (lCoarseNodesPerDim[1] - 1) * lCoarseNodesPerDim[0] + edgeIdx + 1);
528  rowIdx = coordRowIdx * blkSize;
529  coordColumnOffset = ((lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (lFineNodesPerDim[1] - 1) * lFineNodesPerDim[0] + (edgeIdx + 1) * 3);
530  columnOffset = coordColumnOffset * blkSize;
531  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset + edgeLineOffset + (lCoarseNodesPerDim[1] - 2) * faceLineOffset + cornerStencilLength + edgeIdx * edgeStencilLength) * blkSize;
532  for (LO l = 0; l < blkSize; ++l) {
533  for (LO k = 0; k < 3; ++k) {
534  for (LO j = 0; j < 3; ++j) {
535  for (LO i = 0; i < 5; ++i) {
536  entries_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
537  values_h(entryOffset + k * 15 + j * 5 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
538  }
539  }
540  }
541  }
542  for (LO l = 0; l < blkSize; ++l) {
543  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
544  }
545  for (int dim = 0; dim < numDimensions; ++dim) {
546  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
547  }
548  }
549  }
550 
551  // Edges along 1 direction
552  if (lCoarseNodesPerDim[1] - 2 > 0) {
553  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
554  for (LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[1] - 2; ++edgeIdx) {
555  // Edge 0
556  coordRowIdx = (edgeIdx + 1) * lCoarseNodesPerDim[0];
557  rowIdx = coordRowIdx * blkSize;
558  coordColumnOffset = (edgeIdx + 1) * 3 * lFineNodesPerDim[0];
559  columnOffset = coordColumnOffset * blkSize;
560  entryOffset = (edgeLineOffset + edgeIdx * faceLineOffset) * blkSize;
561  for (LO l = 0; l < blkSize; ++l) {
562  for (LO k = 0; k < 3; ++k) {
563  for (LO j = 0; j < 5; ++j) {
564  for (LO i = 0; i < 3; ++i) {
565  entries_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
566  values_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i + 2];
567  }
568  }
569  }
570  }
571  for (LO l = 0; l < blkSize; ++l) {
572  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
573  }
574  for (int dim = 0; dim < numDimensions; ++dim) {
575  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
576  }
577 
578  // Edge 1
579  coordRowIdx = ((edgeIdx + 1) * lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
580  rowIdx = coordRowIdx * blkSize;
581  coordColumnOffset = ((edgeIdx + 1) * 3 * lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
582  columnOffset = coordColumnOffset * blkSize;
583  entryOffset = (edgeLineOffset + edgeIdx * faceLineOffset + edgeStencilLength + (lCoarseNodesPerDim[0] - 2) * faceStencilLength) * blkSize;
584  for (LO l = 0; l < blkSize; ++l) {
585  for (LO k = 0; k < 3; ++k) {
586  for (LO j = 0; j < 5; ++j) {
587  for (LO i = 0; i < 3; ++i) {
588  entries_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
589  values_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i];
590  }
591  }
592  }
593  }
594  for (LO l = 0; l < blkSize; ++l) {
595  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
596  }
597  for (int dim = 0; dim < numDimensions; ++dim) {
598  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
599  }
600 
601  // Edge 2
602  coordRowIdx = ((lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (edgeIdx + 1) * lCoarseNodesPerDim[0]);
603  rowIdx = coordRowIdx * blkSize;
604  coordColumnOffset = ((lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (edgeIdx + 1) * 3 * lFineNodesPerDim[0]);
605  columnOffset = coordColumnOffset * blkSize;
606  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset + edgeLineOffset + edgeIdx * faceLineOffset) * blkSize;
607  for (LO l = 0; l < blkSize; ++l) {
608  for (LO k = 0; k < 3; ++k) {
609  for (LO j = 0; j < 5; ++j) {
610  for (LO i = 0; i < 3; ++i) {
611  entries_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
612  values_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i + 2];
613  }
614  }
615  }
616  }
617  for (LO l = 0; l < blkSize; ++l) {
618  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
619  }
620  for (int dim = 0; dim < numDimensions; ++dim) {
621  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
622  }
623 
624  // Edge 3
625  coordRowIdx = ((lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (edgeIdx + 1) * lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
626  rowIdx = coordRowIdx * blkSize;
627  coordColumnOffset = ((lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (edgeIdx + 1) * 3 * lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
628  columnOffset = coordColumnOffset * blkSize;
629  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset + edgeLineOffset + edgeIdx * faceLineOffset + edgeStencilLength + (lCoarseNodesPerDim[0] - 2) * faceStencilLength) * blkSize;
630  for (LO l = 0; l < blkSize; ++l) {
631  for (LO k = 0; k < 3; ++k) {
632  for (LO j = 0; j < 5; ++j) {
633  for (LO i = 0; i < 3; ++i) {
634  entries_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
635  values_h(entryOffset + k * 15 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
636  }
637  }
638  }
639  }
640  for (LO l = 0; l < blkSize; ++l) {
641  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
642  }
643  for (int dim = 0; dim < numDimensions; ++dim) {
644  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
645  }
646  }
647  }
648 
649  // Edges along 2 direction
650  if (lCoarseNodesPerDim[2] - 2 > 0) {
651  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
652  for (LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[2] - 2; ++edgeIdx) {
653  // Edge 0
654  coordRowIdx = (edgeIdx + 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
655  rowIdx = coordRowIdx * blkSize;
656  coordColumnOffset = (edgeIdx + 1) * 3 * lFineNodesPerDim[1] * lFineNodesPerDim[0];
657  columnOffset = coordColumnOffset * blkSize;
658  entryOffset = (facePlaneOffset + edgeIdx * interiorPlaneOffset) * blkSize;
659  for (LO l = 0; l < blkSize; ++l) {
660  for (LO k = 0; k < 5; ++k) {
661  for (LO j = 0; j < 3; ++j) {
662  for (LO i = 0; i < 3; ++i) {
663  entries_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i) * blkSize + l;
664  values_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i + 2];
665  }
666  }
667  }
668  }
669  for (LO l = 0; l < blkSize; ++l) {
670  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
671  }
672  for (int dim = 0; dim < numDimensions; ++dim) {
673  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
674  }
675 
676  // Edge 1
677  coordRowIdx = ((edgeIdx + 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
678  rowIdx = coordRowIdx * blkSize;
679  coordColumnOffset = ((edgeIdx + 1) * 3 * lFineNodesPerDim[1] * lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
680  columnOffset = coordColumnOffset * blkSize;
681  entryOffset = (facePlaneOffset + faceLineOffset - edgeStencilLength + edgeIdx * interiorPlaneOffset) * blkSize;
682  for (LO l = 0; l < blkSize; ++l) {
683  for (LO k = 0; k < 5; ++k) {
684  for (LO j = 0; j < 3; ++j) {
685  for (LO i = 0; i < 3; ++i) {
686  entries_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i - 2) * blkSize + l;
687  values_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i];
688  }
689  }
690  }
691  }
692  for (LO l = 0; l < blkSize; ++l) {
693  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
694  }
695  for (int dim = 0; dim < numDimensions; ++dim) {
696  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
697  }
698 
699  // Edge 2
700  coordRowIdx = ((edgeIdx + 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (lCoarseNodesPerDim[1] - 1) * lCoarseNodesPerDim[0]);
701  rowIdx = coordRowIdx * blkSize;
702  coordColumnOffset = ((edgeIdx + 1) * 3 * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (lFineNodesPerDim[1] - 1) * lFineNodesPerDim[0]);
703  columnOffset = coordColumnOffset * blkSize;
704  entryOffset = (facePlaneOffset + edgeIdx * interiorPlaneOffset + faceLineOffset + (lCoarseNodesPerDim[1] - 2) * interiorLineOffset) * blkSize;
705  for (LO l = 0; l < blkSize; ++l) {
706  for (LO k = 0; k < 5; ++k) {
707  for (LO j = 0; j < 3; ++j) {
708  for (LO i = 0; i < 3; ++i) {
709  entries_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
710  values_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i + 2];
711  }
712  }
713  }
714  }
715  for (LO l = 0; l < blkSize; ++l) {
716  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
717  }
718  for (int dim = 0; dim < numDimensions; ++dim) {
719  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
720  }
721 
722  // Edge 3
723  coordRowIdx = ((edgeIdx + 2) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] - 1);
724  rowIdx = coordRowIdx * blkSize;
725  coordColumnOffset = ((edgeIdx + 1) * 3 * lFineNodesPerDim[1] * lFineNodesPerDim[0] + lFineNodesPerDim[1] * lFineNodesPerDim[0] - 1);
726  columnOffset = coordColumnOffset * blkSize;
727  entryOffset = (facePlaneOffset + (edgeIdx + 1) * interiorPlaneOffset - edgeStencilLength) * blkSize;
728  for (LO l = 0; l < blkSize; ++l) {
729  for (LO k = 0; k < 5; ++k) {
730  for (LO j = 0; j < 3; ++j) {
731  for (LO i = 0; i < 3; ++i) {
732  entries_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
733  values_h(entryOffset + k * 9 + j * 3 + i + edgeStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
734  }
735  }
736  }
737  }
738  for (LO l = 0; l < blkSize; ++l) {
739  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength * (l + 1);
740  }
741  for (int dim = 0; dim < numDimensions; ++dim) {
742  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
743  }
744  }
745  }
746 
747  // Faces in 0-1 plane
748  if ((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[1] - 2 > 0)) {
749  Array<LO> gridIdx(3);
750  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
751  for (LO faceIdx = 0; faceIdx < (lCoarseNodesPerDim[1] - 2) * (lCoarseNodesPerDim[0] - 2); ++faceIdx) {
752  // Face 0
753  coordRowIdx = ((gridIdx[1] + 1) * lCoarseNodesPerDim[0] + gridIdx[0] + 1);
754  rowIdx = coordRowIdx * blkSize;
755  coordColumnOffset = 3 * ((gridIdx[1] + 1) * lFineNodesPerDim[0] + gridIdx[0] + 1);
756  columnOffset = coordColumnOffset * blkSize;
757  entryOffset = (edgeLineOffset + edgeStencilLength + gridIdx[1] * faceLineOffset + gridIdx[0] * faceStencilLength) * blkSize;
758  for (LO l = 0; l < blkSize; ++l) {
759  for (LO k = 0; k < 3; ++k) {
760  for (LO j = 0; j < 5; ++j) {
761  for (LO i = 0; i < 5; ++i) {
762  entries_h(entryOffset + k * 25 + j * 5 + i + faceStencilLength * l) = columnOffset + (k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
763  values_h(entryOffset + k * 25 + j * 5 + i + faceStencilLength * l) = coeffs[k + 2] * coeffs[j] * coeffs[i];
764  }
765  }
766  }
767  }
768  for (LO l = 0; l < blkSize; ++l) {
769  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
770  }
771  for (int dim = 0; dim < numDimensions; ++dim) {
772  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
773  }
774 
775  // Face 1
776  coordRowIdx += (lCoarseNodesPerDim[2] - 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0];
777  rowIdx = coordRowIdx * blkSize;
778  coordColumnOffset += (lFineNodesPerDim[2] - 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0];
779  columnOffset = coordColumnOffset * blkSize;
780  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2) * interiorPlaneOffset) * blkSize;
781  for (LO l = 0; l < blkSize; ++l) {
782  for (LO k = 0; k < 3; ++k) {
783  for (LO j = 0; j < 5; ++j) {
784  for (LO i = 0; i < 5; ++i) {
785  entries_h(entryOffset + k * 25 + j * 5 + i + faceStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
786  values_h(entryOffset + k * 25 + j * 5 + i + faceStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
787  }
788  }
789  }
790  }
791  for (LO l = 0; l < blkSize; ++l) {
792  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
793  }
794  for (int dim = 0; dim < numDimensions; ++dim) {
795  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
796  }
797 
798  // Last step in the loop
799  // update the grid indices
800  // for next grid point
801  ++gridIdx[0];
802  if (gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
803  gridIdx[0] = 0;
804  ++gridIdx[1];
805  }
806  }
807  }
808 
809  // Faces in 0-2 plane
810  if ((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
811  Array<LO> gridIdx(3);
812  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
813  for (LO faceIdx = 0; faceIdx < (lCoarseNodesPerDim[2] - 2) * (lCoarseNodesPerDim[0] - 2); ++faceIdx) {
814  // Face 0
815  coordRowIdx = ((gridIdx[2] + 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (gridIdx[0] + 1));
816  rowIdx = coordRowIdx * blkSize;
817  coordColumnOffset = ((gridIdx[2] + 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + gridIdx[0] + 1) * 3;
818  columnOffset = coordColumnOffset * blkSize;
819  entryOffset = (facePlaneOffset + gridIdx[2] * interiorPlaneOffset + edgeStencilLength + gridIdx[0] * faceStencilLength) * blkSize;
820  for (LO l = 0; l < blkSize; ++l) {
821  for (LO k = 0; k < 5; ++k) {
822  for (LO j = 0; j < 3; ++j) {
823  for (LO i = 0; i < 5; ++i) {
824  entries_h(entryOffset + k * 15 + j * 5 + i + faceStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i - 2) * blkSize + l;
825  values_h(entryOffset + k * 15 + j * 5 + i + faceStencilLength * l) = coeffs[k] * coeffs[j + 2] * coeffs[i];
826  }
827  }
828  }
829  }
830  for (LO l = 0; l < blkSize; ++l) {
831  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
832  }
833  for (int dim = 0; dim < numDimensions; ++dim) {
834  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
835  }
836 
837  // Face 1
838  coordRowIdx += (lCoarseNodesPerDim[1] - 1) * lCoarseNodesPerDim[0];
839  rowIdx = coordRowIdx * blkSize;
840  coordColumnOffset += (lFineNodesPerDim[1] - 1) * lFineNodesPerDim[0];
841  columnOffset = coordColumnOffset * blkSize;
842  entryOffset += (faceLineOffset + (lCoarseNodesPerDim[1] - 2) * interiorLineOffset) * blkSize;
843  for (LO l = 0; l < blkSize; ++l) {
844  for (LO k = 0; k < 5; ++k) {
845  for (LO j = 0; j < 3; ++j) {
846  for (LO i = 0; i < 5; ++i) {
847  entries_h(entryOffset + k * 15 + j * 5 + i + faceStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
848  values_h(entryOffset + k * 15 + j * 5 + i + faceStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
849  }
850  }
851  }
852  }
853  for (LO l = 0; l < blkSize; ++l) {
854  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
855  }
856  for (int dim = 0; dim < numDimensions; ++dim) {
857  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
858  }
859 
860  // Last step in the loop
861  // update the grid indices
862  // for next grid point
863  ++gridIdx[0];
864  if (gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
865  gridIdx[0] = 0;
866  ++gridIdx[2];
867  }
868  }
869  }
870 
871  // Faces in 1-2 plane
872  if ((lCoarseNodesPerDim[1] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
873  Array<LO> gridIdx(3);
874  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
875  for (LO faceIdx = 0; faceIdx < (lCoarseNodesPerDim[2] - 2) * (lCoarseNodesPerDim[1] - 2); ++faceIdx) {
876  // Face 0
877  coordRowIdx = ((gridIdx[2] + 1) * lCoarseNodesPerDim[1] * lCoarseNodesPerDim[0] + (gridIdx[1] + 1) * lCoarseNodesPerDim[0]);
878  rowIdx = coordRowIdx * blkSize;
879  coordColumnOffset = ((gridIdx[2] + 1) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (gridIdx[1] + 1) * lFineNodesPerDim[0]) * 3;
880  columnOffset = coordColumnOffset * blkSize;
881  entryOffset = (facePlaneOffset + gridIdx[2] * interiorPlaneOffset + faceLineOffset + gridIdx[1] * interiorLineOffset) * blkSize;
882  for (LO l = 0; l < blkSize; ++l) {
883  for (LO k = 0; k < 5; ++k) {
884  for (LO j = 0; j < 5; ++j) {
885  for (LO i = 0; i < 3; ++i) {
886  entries_h(entryOffset + k * 15 + j * 3 + i + faceStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i) * blkSize + l;
887  values_h(entryOffset + k * 15 + j * 3 + i + faceStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i + 2];
888  }
889  }
890  }
891  }
892  for (LO l = 0; l < blkSize; ++l) {
893  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
894  }
895  for (int dim = 0; dim < numDimensions; ++dim) {
896  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
897  }
898 
899  // Face 1
900  coordRowIdx += (lCoarseNodesPerDim[0] - 1);
901  rowIdx = coordRowIdx * blkSize;
902  coordColumnOffset += (lFineNodesPerDim[0] - 1);
903  columnOffset = coordColumnOffset * blkSize;
904  entryOffset += (faceStencilLength + (lCoarseNodesPerDim[0] - 2) * interiorStencilLength) * blkSize;
905  for (LO l = 0; l < blkSize; ++l) {
906  for (LO k = 0; k < 5; ++k) {
907  for (LO j = 0; j < 5; ++j) {
908  for (LO i = 0; i < 3; ++i) {
909  entries_h(entryOffset + k * 15 + j * 3 + i + faceStencilLength * l) = columnOffset + ((k - 2) * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (j - 2) * lFineNodesPerDim[0] + i - 2) * blkSize + l;
910  values_h(entryOffset + k * 15 + j * 3 + i + faceStencilLength * l) = coeffs[k] * coeffs[j] * coeffs[i];
911  }
912  }
913  }
914  }
915  for (LO l = 0; l < blkSize; ++l) {
916  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength * (l + 1);
917  }
918  for (int dim = 0; dim < numDimensions; ++dim) {
919  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
920  }
921 
922  // Last step in the loop
923  // update the grid indices
924  // for next grid point
925  ++gridIdx[1];
926  if (gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
927  gridIdx[1] = 0;
928  ++gridIdx[2];
929  }
930  }
931  }
932 
933  if (numInteriors > 0) {
934  // Allocate and compute arrays
935  // containing column offsets
936  // and values associated with
937  // interior points
938  LO countRowEntries = 0;
939  Array<LO> coordColumnOffsets(125);
940  for (LO k = -2; k < 3; ++k) {
941  for (LO j = -2; j < 3; ++j) {
942  for (LO i = -2; i < 3; ++i) {
943  coordColumnOffsets[countRowEntries] = k * lFineNodesPerDim[1] * lFineNodesPerDim[0] + j * lFineNodesPerDim[0] + i;
944  ++countRowEntries;
945  }
946  }
947  }
948 
949  LO countValues = 0;
950  Array<SC> interiorValues(125);
951  for (LO k = 0; k < 5; ++k) {
952  for (LO j = 0; j < 5; ++j) {
953  for (LO i = 0; i < 5; ++i) {
954  interiorValues[countValues] = coeffs[k] * coeffs[j] * coeffs[i];
955  ++countValues;
956  }
957  }
958  }
959 
960  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
961  Array<LO> gridIdx(3);
962  for (LO interiorIdx = 0; interiorIdx < numInteriors; ++interiorIdx) {
963  coordRowIdx = ((gridIdx[2] + 1) * lCoarseNodesPerDim[0] * lCoarseNodesPerDim[1] + (gridIdx[1] + 1) * lCoarseNodesPerDim[0] + gridIdx[0] + 1);
964  rowIdx = coordRowIdx * blkSize;
965  coordColumnOffset = ((gridIdx[2] + 1) * 3 * lFineNodesPerDim[1] * lFineNodesPerDim[0] + (gridIdx[1] + 1) * 3 * lFineNodesPerDim[0] + (gridIdx[0] + 1) * 3);
966  columnOffset = coordColumnOffset * blkSize;
967  entryOffset = (facePlaneOffset + faceLineOffset + faceStencilLength + gridIdx[2] * interiorPlaneOffset + gridIdx[1] * interiorLineOffset + gridIdx[0] * interiorStencilLength) * blkSize;
968  for (LO l = 0; l < blkSize; ++l) {
969  row_map_h(rowIdx + 1 + l) = entryOffset + interiorStencilLength * (l + 1);
970  }
971  // Fill the column indices
972  // and values in the approproate
973  // views.
974  for (LO l = 0; l < blkSize; ++l) {
975  for (LO entryIdx = 0; entryIdx < interiorStencilLength; ++entryIdx) {
976  entries_h(entryOffset + entryIdx + interiorStencilLength * l) = columnOffset + coordColumnOffsets[entryIdx] * blkSize + l;
977  values_h(entryOffset + entryIdx + interiorStencilLength * l) = interiorValues[entryIdx];
978  }
979  }
980  for (int dim = 0; dim < numDimensions; ++dim) {
981  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
982  }
983 
984  // Last step in the loop
985  // update the grid indices
986  // for next grid point
987  ++gridIdx[0];
988  if (gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
989  gridIdx[0] = 0;
990  ++gridIdx[1];
991  if (gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
992  gridIdx[1] = 0;
993  ++gridIdx[2];
994  }
995  }
996  }
997  }
998 
999  Kokkos::deep_copy(row_map, row_map_h);
1000  Kokkos::deep_copy(entries, entries_h);
1001  Kokkos::deep_copy(values, values_h);
1002 
1003  local_graph_type localGraph(entries, row_map);
1004  local_matrix_type localR("R", numCols, values, localGraph);
1005 
1006  R = MatrixFactory::Build(localR, // the local data
1007  rowMap, // rowMap
1008  A->getRowMap(), // colMap
1009  A->getRowMap(), // domainMap == colMap
1010  rowMap, // rangeMap == rowMap
1011  Teuchos::null); // params for optimized construction
1012 
1013 } // Build3D()
1014 
1015 } // namespace MueLu
1016 
1017 #define MUELU_REGIONRFACTORY_SHORT
1018 #endif // MUELU_REGIONRFACTORY_DEF_HPP
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:13
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
Teuchos::basic_FancyOStream::setShowProcRank
basic_FancyOStream & setShowProcRank(const bool showProcRank)
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
Teuchos::basic_oblackholestream
Xpetra_CrsGraphFactory.hpp
LO
LocalOrdinal LO
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Xpetra_Matrix.hpp
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:15
MueLu_RegionRFactory_decl.hpp
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:34
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:14
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
Teuchos::ParameterList::isSublist
bool isSublist(const std::string &name) const
MueLu::RegionRFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_RegionRFactory_def.hpp:56
MueLu::Utilities::Transpose
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 > &params=Teuchos::null)
Definition: MueLu_Utilities_def.hpp:496
Teuchos::basic_FancyOStream::setShowAllFrontMatter
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:120
MueLu::RegionRFactory::Build3D
void Build3D(const int numDimensions, Array< LO > &lFineNodesPerDim, const RCP< Matrix > &A, const RCP< realvaluedmultivector_type > &fineCoordinates, RCP< Matrix > &R, RCP< realvaluedmultivector_type > &coarseCoordinates, Array< LO > &lCoarseNodesPerDim) const
Definition: MueLu_RegionRFactory_def.hpp:156
Teuchos::ParameterList::sublist
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
MueLu::RegionRFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_RegionRFactory_def.hpp:45
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:51
MueLu::RegionRFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Input.
Definition: MueLu_RegionRFactory_def.hpp:25
Teuchos::toString
std::string toString(const T &t)