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