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