MueLu Version of the Day
Loading...
Searching...
No Matches
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
21namespace MueLu {
22
23template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
24RCP<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
43template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
44void 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
54template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
55void 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
154template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
155void 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_device_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::host_mirror_type row_map_h = Kokkos::create_mirror_view(row_map);
244
245 entries_type entries(Kokkos::ViewAllocateWithoutInitializing("entries"), nnz);
246 typename entries_type::host_mirror_type entries_h = Kokkos::create_mirror_view(entries);
247
248 values_type values(Kokkos::ViewAllocateWithoutInitializing("values"), nnz);
249 typename values_type::host_mirror_type 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
MueLu_RegionRFactory_decl.hpp
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition MueLu_UseDefaultTypes.hpp:13
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition MueLu_UseDefaultTypes.hpp:14
Node
MueLu::DefaultNode Node
Definition MueLu_UseDefaultTypes.hpp:15
MueLu::Level
Class that holds all level-specific information.
Definition MueLu_Level.hpp:63
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition MueLu_Level.cpp:51
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition MueLu_Level.hpp:120
MueLu::NoFactory::get
static const NoFactory * get()
Definition MueLu_NoFactory.cpp:34
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
MueLu::RegionRFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition MueLu_RegionRFactory_def.hpp:45
MueLu::RegionRFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition MueLu_RegionRFactory_def.hpp:56
MueLu::RegionRFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Input.
Definition MueLu_RegionRFactory_def.hpp:25
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:178
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42