MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_GeoInterpFactory_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_GEOINTERPFACTORY_DEF_HPP
11#define MUELU_GEOINTERPFACTORY_DEF_HPP
12
13#include <iostream>
14#include <cmath>
15
16#include <Teuchos_SerialDenseMatrix.hpp>
17
18#include <Xpetra_Map.hpp>
19#include <Xpetra_MapFactory.hpp>
20#include <Xpetra_Matrix.hpp>
21#include <Xpetra_MultiVector.hpp>
22#include <Xpetra_MultiVectorFactory.hpp>
23#include <Xpetra_VectorFactory.hpp>
24#include <Xpetra_CrsMatrixWrap.hpp>
25
26#include "MueLu_GeoInterpFactory_decl.hpp"
27#include <MueLu_Level.hpp>
28
29namespace MueLu {
30
31template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
32GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GeoInterpFactory() {
33 GetOStream(Runtime1) << "I constructed a GeoInterpFactory object... Nothing else to do here." << std::endl;
34}
35
36template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
37GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~GeoInterpFactory() {
38 // Should be empty. All destruction should be handled by Level-based get stuff and RCP
39}
40
41template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
42void GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::DeclareInput(Level &fineLevel, Level &coarseLevel) const {
43 Input(fineLevel, "A");
44 Input(fineLevel, "A00");
45 Input(fineLevel, "A10");
46 Input(fineLevel, "A20");
47
48 Input(fineLevel, "VElementList");
49 Input(fineLevel, "PElementList");
50 Input(fineLevel, "MElementList");
51
52 Input(coarseLevel, "VElementList");
53 Input(coarseLevel, "PElementList");
54 Input(coarseLevel, "MElementList");
55
56 /*
57 coarseLevel.DeclareInput("VElementList",coarseLevel.GetFactoryManager()->GetFactory("VElementList").get(),this);
58 coarseLevel.DeclareInput("PElementList",coarseLevel.GetFactoryManager()->GetFactory("PElementList").get(),this);
59 coarseLevel.DeclareInput("MElementList",coarseLevel.GetFactoryManager()->GetFactory("PElementList").get(),this);
60
61 fineLevel.DeclareInput("VElementList",fineLevel.GetFactoryManager()->GetFactory("VElementList").get(),this);
62 fineLevel.DeclareInput("PElementList",fineLevel.GetFactoryManager()->GetFactory("PElementList").get(),this);
63 fineLevel.DeclareInput("MElementList",fineLevel.GetFactoryManager()->GetFactory("PElementList").get(),this);
64 */
65
66 // currentLevel.DeclareInput(varName_,factory_,this);
67}
68
69template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70void GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Level &fineLevel, Level &coarseLevel) const {
71 GetOStream(Runtime1) << "Starting 'build' routine...\n";
72
73 // This will create a list of elements on the coarse grid with a
74 // predictable structure, as well as modify the fine grid list of
75 // elements, if necessary (i.e. if fineLevel.GetLevelID()==0);
76 // BuildCoarseGrid(fineLevel,coarseLevel);
77
78 // This will actually build our prolongator P
79 return BuildP(fineLevel, coarseLevel);
80}
81
82template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
83void GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level &fineLevel, Level &coarseLevel) const {
84 typedef Teuchos::SerialDenseMatrix<GO, GO> SerialDenseMatrixType;
85
86 GetOStream(Runtime1) << "Starting 'BuildP' routine...\n";
87
88 // DEBUG
89 // Teuchos::FancyOStream fout(*GetOStream(Runtime1));
90 // fineLevel.print(fout,Teuchos::VERB_HIGH);
91
92 // Get finegrid element lists
93 RCP<SerialDenseMatrixType> fineElementPDOFs = Get<RCP<SerialDenseMatrixType> >(fineLevel, "PElementList");
94 RCP<SerialDenseMatrixType> fineElementVDOFs = Get<RCP<SerialDenseMatrixType> >(fineLevel, "VElementList");
95 RCP<SerialDenseMatrixType> fineElementMDOFs = Get<RCP<SerialDenseMatrixType> >(fineLevel, "MElementList");
96
97 // DEBUG
98 GetOStream(Runtime1) << "done getting fine level elements...\n";
99 GetOStream(Runtime1) << "getting coarse level elements...\n";
100 // coarseLevel.print(fout,Teuchos::VERB_HIGH);
101
102 // Get coarse grid element lists
103 RCP<Teuchos::SerialDenseMatrix<GO, GO> > coarseElementVDOFs,
104 coarseElementPDOFs,
105 coarseElementMDOFs;
106
107 coarseLevel.Get("VElementList", coarseElementVDOFs, coarseLevel.GetFactoryManager()->GetFactory("VElementList").get());
108 coarseLevel.Get("PElementList", coarseElementPDOFs, coarseLevel.GetFactoryManager()->GetFactory("PElementList").get());
109 coarseLevel.Get("MElementList", coarseElementMDOFs, coarseLevel.GetFactoryManager()->GetFactory("MElementList").get());
110
111 GetOStream(Runtime1) << "computing various numbers...\n";
112 // Number of elements?
113 GO totalFineElements = fineElementMDOFs->numRows();
114 LO nFineElements = (int)sqrt(totalFineElements);
115 GO totalCoarseElements = coarseElementMDOFs->numRows();
116 LO nCoarseElements = (int)sqrt(totalCoarseElements);
117
118 // Set sizes for *COARSE GRID*
119 GO nM = (2 * nCoarseElements + 1) * (2 * nCoarseElements + 1);
120 GO nV = 2 * nM;
121 GO nP = (nCoarseElements + 1) * (nCoarseElements + 1);
122
123 // Get the row maps for the Ps
124 RCP<Matrix> fineA00 = Get<RCP<Matrix> >(fineLevel, "A00");
125 RCP<Matrix> fineA10 = Get<RCP<Matrix> >(fineLevel, "A10");
126 RCP<Matrix> fineA20 = Get<RCP<Matrix> >(fineLevel, "A20");
127
128 GetOStream(Runtime1) << "creating coarse grid maps...\n";
129
130 RCP<const Map> rowMapforPV = fineA00->getRowMap();
131 RCP<const Map> rowMapforPP = fineA10->getRowMap();
132 RCP<const Map> rowMapforPM = fineA20->getRowMap();
133
134 GO fNV = rowMapforPV->getGlobalNumElements();
135 GO fNP = rowMapforPP->getGlobalNumElements();
136 GO fNM = rowMapforPM->getGlobalNumElements();
137
138 // Get the comm for the maps
139 RCP<const Teuchos::Comm<int> > comm = rowMapforPV->getComm();
140
141 // Create rowMap for P
142 RCP<Matrix> FineA = Factory::Get<RCP<Matrix> >(fineLevel, "A");
143 RCP<const Map> rowMapforP = FineA->getRowMap();
144
145 // Create colMaps for the coarse grid
146 RCP<const Map> colMapforPV = Xpetra::MapFactory<LO, GO>::createUniformContigMap(Xpetra::UseTpetra, nV, comm);
147 RCP<const Map> colMapforPP = Xpetra::MapFactory<LO, GO>::createUniformContigMap(Xpetra::UseTpetra, nP, comm);
148 RCP<const Map> colMapforPM = Xpetra::MapFactory<LO, GO>::createUniformContigMap(Xpetra::UseTpetra, nM, comm);
149
150 GetOStream(Runtime1) << "creating coarse grid matrices...\n";
151 // Create our final output Ps for the coarseGrid
152 size_t maxEntriesPerRowV = 9, // No overlap of VX and VY
153 maxEntriesPerRowP = 4,
154 maxEntriesPerRowM = 9;
155
156 RCP<Matrix> P = rcp(new CrsMatrixWrap(rowMapforP, maxEntriesPerRowV));
157 RCP<Matrix> PV = rcp(new CrsMatrixWrap(rowMapforPV, maxEntriesPerRowV));
158 RCP<Matrix> PP = rcp(new CrsMatrixWrap(rowMapforPP, maxEntriesPerRowP));
159 RCP<Matrix> PM = rcp(new CrsMatrixWrap(rowMapforPM, maxEntriesPerRowM));
160
161 //*****************************************************************/
162 //
163 // All 25 fine grid dofs are completely determined by the coarse
164 // grid element in which they reside! So I should loop over coarse
165 // grid elements and build 25 rows at a time! If that's not
166 // ridiculous... I just have to be careful about duplicates on
167 // future elements! But duplicates are easy - just the bottom and
168 // left edges.
169 //
170 //
171 // Looking at a fine grid patch, define the following Local-Global
172 // relationship (magnetics as an example):
173 //
174 // Bottom-Left Corner:
175 // 0 -> (*fineElementMDOFs)(fineElement[0],0)
176 //
177 // Bottom Edge:
178 // 1 -> (*fineElementMDOFs)(fineElement[0],4)
179 // 2 -> (*fineElementMDOFs)(fineElement[0],2)
180 // 3 -> (*fineElementMDOFs)(fineElement[1],4)
181 // 4 -> (*fineElementMDOFs)(fineElement[1],2)
182 //
183 // Left Edge:
184 // 5 -> (*fineElementMDOFs)(fineElement[0],7)
185 // 6 -> (*fineElementMDOFs)(fineElement[0],3)
186 // 7 -> (*fineElementMDOFs)(fineElement[2],7)
187 // 8 -> (*fineElementMDOFs)(fineElement[2],3)
188 //
189 // All the rest:
190 // 9 -> (*fineElementMDOFs)(fineElement[3],0)
191 // 10 -> (*fineElementMDOFs)(fineElement[3],1)
192 // 11 -> (*fineElementMDOFs)(fineElement[3],2)
193 // 12 -> (*fineElementMDOFs)(fineElement[3],3)
194 // 13 -> (*fineElementMDOFs)(fineElement[0],5)
195 // 14 -> (*fineElementMDOFs)(fineElement[0],6)
196 // 15 -> (*fineElementMDOFs)(fineElement[1],5)
197 // 16 -> (*fineElementMDOFs)(fineElement[1],6)
198 // 17 -> (*fineElementMDOFs)(fineElement[2],5)
199 // 18 -> (*fineElementMDOFs)(fineElement[2],6)
200 // 19 -> (*fineElementMDOFs)(fineElement[3],5)
201 // 20 -> (*fineElementMDOFs)(fineElement[3],6)
202 // 21 -> (*fineElementMDOFs)(fineElement[0],8)
203 // 22 -> (*fineElementMDOFs)(fineElement[1],8)
204 // 23 -> (*fineElementMDOFs)(fineElement[2],8)
205 // 24 -> (*fineElementMDOFs)(fineElement[3],8)
206 //
207 //*****************************************************************/
208
209 size_t nnz = 0; // Just to make my copy-paste life easier...
210 Teuchos::ArrayRCP<GO> colPtrV(maxEntriesPerRowV, 0);
211 Teuchos::ArrayRCP<GO> colPtrM(maxEntriesPerRowM, 0);
212 Teuchos::ArrayRCP<SC> valPtrM(maxEntriesPerRowM, 0.);
213
214 Teuchos::ArrayRCP<GO> colPtrP(maxEntriesPerRowP, 0);
215 Teuchos::ArrayRCP<SC> valPtrP(maxEntriesPerRowP, 0.);
216
217 // About which fine-grid elements do we care?
218 GO fineElement[4] = {0, 1, nFineElements, nFineElements + 1};
219
220 GetOStream(Runtime1) << "start building matrices...\n";
221 GetOStream(Runtime1) << "nCoarseElements = " << nCoarseElements << std::endl;
222
223 for (GO coarseElement = 0; coarseElement < totalCoarseElements; coarseElement++) {
224 // We don't really care what is shared with future elements -
225 // we know a priori what needs to be filled for the element
226 // we're dealing with, depending of if it it's on the bottom
227 // edge, the left edge, or in the middle. Thoses should be the
228 // only cases we care about, and they should require a
229 // superset of the work required for an interior node.
230
231 // if (CoarseElement is on bottom edge)
232 if (coarseElement < nCoarseElements) {
233 // fill in the bottom edge of the element patch
234 // FP = 1
235 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
236 valPtrM[0] = 0.375;
237 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
238 valPtrM[1] = -0.125;
239 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 4);
240 valPtrM[2] = 0.75;
241
242 nnz = 3;
243 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 4), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
244
245 colPtrV[0] = 2 * colPtrM[0];
246 colPtrV[1] = 2 * colPtrM[1];
247 colPtrV[2] = 2 * colPtrM[2];
248
249 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 8), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
250
251 colPtrV[0] = 2 * colPtrM[0] + 1;
252 colPtrV[1] = 2 * colPtrM[1] + 1;
253 colPtrV[2] = 2 * colPtrM[2] + 1;
254
255 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 9), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
256
257 // FPr = 1
258 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 0);
259 valPtrP[0] = 0.5;
260 colPtrP[1] = (*coarseElementPDOFs)(coarseElement, 1);
261 valPtrP[1] = 0.5;
262
263 nnz = 2;
264 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[0], 1), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
265
266 // FP = 2
267 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 4);
268 valPtrM[0] = 1.0;
269
270 nnz = 1;
271 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 1), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
272
273 colPtrV[0] = 2 * colPtrM[0];
274
275 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 2), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
276
277 colPtrV[0] = 2 * colPtrM[0] + 1;
278
279 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 3), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
280
281 // FPr = 2
282 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 1);
283 valPtrP[0] = 1.0;
284
285 nnz = 1;
286 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[1], 1), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
287
288 // FP = 3
289 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
290 valPtrM[0] = -0.125;
291 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
292 valPtrM[1] = 0.375;
293 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 4);
294 valPtrM[2] = 0.75;
295
296 nnz = 3;
297 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[1], 4), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
298
299 colPtrV[0] = 2 * colPtrM[0];
300 colPtrV[1] = 2 * colPtrM[1];
301 colPtrV[2] = 2 * colPtrM[2];
302
303 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 8), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
304
305 colPtrV[0] = 2 * colPtrM[0] + 1;
306 colPtrV[1] = 2 * colPtrM[1] + 1;
307 colPtrV[2] = 2 * colPtrM[2] + 1;
308
309 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 9), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
310
311 // FP = 4
312 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 1);
313 valPtrM[0] = 1.0;
314
315 nnz = 1;
316 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[1], 1), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
317
318 colPtrV[0] = 2 * colPtrM[0];
319
320 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 2), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
321
322 colPtrV[0] = 2 * colPtrM[0] + 1;
323
324 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 3), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
325
326 // if (CoarseElement is on the bottom left corner)
327 if (coarseElement == 0) {
328 // fill in the bottom left corner
329 // FP = 0
330 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
331 valPtrM[0] = 1.0;
332
333 nnz = 1;
334 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 0), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
335
336 colPtrV[0] = 2 * colPtrM[0];
337
338 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 0), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
339
340 colPtrV[0] = 2 * colPtrM[0] + 1;
341
342 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 1), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
343
344 // FPr = 0
345 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 0);
346 valPtrP[0] = 1.0;
347
348 nnz = 1;
349 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[0], 0), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
350
351 } // if (coarseElement is on the bottom left corner)
352 } // if (coarseElement is on the bottom edge)
353
354 // if (CoarseElement is on left edge)
355 if (coarseElement % (nCoarseElements) == 0) {
356 // fill in the left edge of the element patch
357 // FP = 5
358 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
359 valPtrM[0] = 0.375;
360 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 3);
361 valPtrM[1] = -0.125;
362 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 7);
363 valPtrM[2] = 0.75;
364
365 nnz = 3;
366 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 7), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
367
368 colPtrV[0] = 2 * colPtrM[0];
369 colPtrV[1] = 2 * colPtrM[1];
370 colPtrV[2] = 2 * colPtrM[2];
371
372 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 14), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
373
374 colPtrV[0] = 2 * colPtrM[0] + 1;
375 colPtrV[1] = 2 * colPtrM[1] + 1;
376 colPtrV[2] = 2 * colPtrM[2] + 1;
377
378 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 15), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
379
380 // FP = 6
381 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 7);
382 valPtrM[0] = 1.0;
383
384 nnz = 1;
385 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 3), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
386
387 colPtrV[0] = 2 * colPtrM[0];
388
389 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 6), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
390
391 colPtrV[0] = 2 * colPtrM[0] + 1;
392
393 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 7), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
394
395 // FP = 7
396 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
397 valPtrM[0] = -0.125;
398 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 3);
399 valPtrM[1] = 0.375;
400 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 7);
401 valPtrM[2] = 0.75;
402
403 nnz = 3;
404 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[2], 7), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
405
406 colPtrV[0] = 2 * colPtrM[0];
407 colPtrV[1] = 2 * colPtrM[1];
408 colPtrV[2] = 2 * colPtrM[2];
409
410 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 14), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
411
412 colPtrV[0] = 2 * colPtrM[0] + 1;
413 colPtrV[1] = 2 * colPtrM[1] + 1;
414 colPtrV[2] = 2 * colPtrM[2] + 1;
415
416 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 15), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
417
418 // FP = 8
419 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 3);
420 valPtrM[0] = 1.0;
421
422 nnz = 1;
423 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[2], 3), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
424
425 colPtrV[0] = 2 * colPtrM[0];
426
427 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 6), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
428
429 colPtrV[0] = 2 * colPtrM[0] + 1;
430
431 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 7), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
432
433 // FPr = 3
434 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 0);
435 valPtrP[0] = 0.5;
436 colPtrP[1] = (*coarseElementPDOFs)(coarseElement, 3);
437 valPtrP[1] = 0.5;
438
439 nnz = 2;
440 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[0], 3), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
441
442 // FPr = 4
443 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 3);
444 valPtrP[0] = 1.0;
445
446 nnz = 1;
447 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[2], 3), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
448
449 } // endif (coarseElement is on left edge)
450
451 // fill in the rest of the patch
452 // FP = 9
453 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 8);
454 valPtrM[0] = 1.0;
455
456 nnz = 1;
457 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 0), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
458
459 colPtrV[0] = 2 * colPtrM[0];
460
461 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 0), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
462
463 colPtrV[0] = 2 * colPtrM[0] + 1;
464
465 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 1), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
466
467 // FP = 10
468 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 5);
469 valPtrM[0] = 1.0;
470
471 nnz = 1;
472 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 1), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
473
474 colPtrV[0] = 2 * colPtrM[0];
475
476 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 2), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
477
478 colPtrV[0] = 2 * colPtrM[0] + 1;
479
480 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 3), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
481
482 // FP = 11
483 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 2);
484 valPtrM[0] = 1.0;
485
486 nnz = 1;
487 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 2), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
488
489 colPtrV[0] = 2 * colPtrM[0];
490
491 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 4), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
492
493 colPtrV[0] = 2 * colPtrM[0] + 1;
494
495 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 5), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
496
497 // FP = 12
498 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 6);
499 valPtrM[0] = 1.0;
500
501 nnz = 1;
502 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 3), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
503
504 colPtrV[0] = 2 * colPtrM[0];
505
506 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 6), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
507
508 colPtrV[0] = 2 * colPtrM[0] + 1;
509
510 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 7), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
511
512 // FP = 13
513 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 4);
514 valPtrM[0] = 0.375;
515 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 6);
516 valPtrM[1] = -0.125;
517 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 8);
518 valPtrM[2] = 0.75;
519
520 nnz = 3;
521 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 5), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
522
523 colPtrV[0] = 2 * colPtrM[0];
524 colPtrV[1] = 2 * colPtrM[1];
525 colPtrV[2] = 2 * colPtrM[2];
526
527 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 10), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
528
529 colPtrV[0] = 2 * colPtrM[0] + 1;
530 colPtrV[1] = 2 * colPtrM[1] + 1;
531 colPtrV[2] = 2 * colPtrM[2] + 1;
532
533 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 11), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
534
535 // FP = 14
536 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 5);
537 valPtrM[0] = -0.125;
538 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 7);
539 valPtrM[1] = 0.375;
540 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 8);
541 valPtrM[2] = 0.75;
542
543 nnz = 3;
544 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 6), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
545
546 colPtrV[0] = 2 * colPtrM[0];
547 colPtrV[1] = 2 * colPtrM[1];
548 colPtrV[2] = 2 * colPtrM[2];
549
550 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 12), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
551
552 colPtrV[0] = 2 * colPtrM[0] + 1;
553 colPtrV[1] = 2 * colPtrM[1] + 1;
554 colPtrV[2] = 2 * colPtrM[2] + 1;
555
556 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 13), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
557
558 // FP = 15
559 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 1);
560 valPtrM[0] = 0.375;
561 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 2);
562 valPtrM[1] = -0.125;
563 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 5);
564 valPtrM[2] = 0.75;
565
566 nnz = 3;
567 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[1], 5), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
568
569 colPtrV[0] = 2 * colPtrM[0];
570 colPtrV[1] = 2 * colPtrM[1];
571 colPtrV[2] = 2 * colPtrM[2];
572
573 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 10), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
574
575 colPtrV[0] = 2 * colPtrM[0] + 1;
576 colPtrV[1] = 2 * colPtrM[1] + 1;
577 colPtrV[2] = 2 * colPtrM[2] + 1;
578
579 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 11), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
580
581 // FP = 16
582 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 5);
583 valPtrM[0] = 0.375;
584 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 7);
585 valPtrM[1] = -0.125;
586 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 8);
587 valPtrM[2] = 0.75;
588
589 nnz = 3;
590 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[1], 6), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
591
592 colPtrV[0] = 2 * colPtrM[0];
593 colPtrV[1] = 2 * colPtrM[1];
594 colPtrV[2] = 2 * colPtrM[2];
595
596 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 12), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
597
598 colPtrV[0] = 2 * colPtrM[0] + 1;
599 colPtrV[1] = 2 * colPtrM[1] + 1;
600 colPtrV[2] = 2 * colPtrM[2] + 1;
601
602 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 13), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
603
604 // FP = 17
605 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 4);
606 valPtrM[0] = -0.125;
607 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 6);
608 valPtrM[1] = 0.375;
609 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 8);
610 valPtrM[2] = 0.75;
611
612 nnz = 3;
613 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[2], 5), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
614
615 colPtrV[0] = 2 * colPtrM[0];
616 colPtrV[1] = 2 * colPtrM[1];
617 colPtrV[2] = 2 * colPtrM[2];
618
619 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 10), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
620
621 colPtrV[0] = 2 * colPtrM[0] + 1;
622 colPtrV[1] = 2 * colPtrM[1] + 1;
623 colPtrV[2] = 2 * colPtrM[2] + 1;
624
625 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 11), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
626
627 // FP = 18
628 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 2);
629 valPtrM[0] = -0.125;
630 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 3);
631 valPtrM[1] = 0.375;
632 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 6);
633 valPtrM[2] = 0.75;
634
635 nnz = 3;
636 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[2], 6), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
637
638 colPtrV[0] = 2 * colPtrM[0];
639 colPtrV[1] = 2 * colPtrM[1];
640 colPtrV[2] = 2 * colPtrM[2];
641
642 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 12), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
643
644 colPtrV[0] = 2 * colPtrM[0] + 1;
645 colPtrV[1] = 2 * colPtrM[1] + 1;
646 colPtrV[2] = 2 * colPtrM[2] + 1;
647
648 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 13), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
649
650 // FP = 19
651 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 1);
652 valPtrM[0] = -0.125;
653 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 2);
654 valPtrM[1] = 0.375;
655 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 5);
656 valPtrM[2] = 0.75;
657
658 nnz = 3;
659 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 5), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
660
661 colPtrV[0] = 2 * colPtrM[0];
662 colPtrV[1] = 2 * colPtrM[1];
663 colPtrV[2] = 2 * colPtrM[2];
664
665 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 10), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
666
667 colPtrV[0] = 2 * colPtrM[0] + 1;
668 colPtrV[1] = 2 * colPtrM[1] + 1;
669 colPtrV[2] = 2 * colPtrM[2] + 1;
670
671 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 11), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
672
673 // FP = 20
674 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 2);
675 valPtrM[0] = 0.375;
676 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 3);
677 valPtrM[1] = -0.125;
678 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 6);
679 valPtrM[2] = 0.75;
680
681 nnz = 3;
682 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 6), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
683
684 colPtrV[0] = 2 * colPtrM[0];
685 colPtrV[1] = 2 * colPtrM[1];
686 colPtrV[2] = 2 * colPtrM[2];
687
688 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 12), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
689
690 colPtrV[0] = 2 * colPtrM[0] + 1;
691 colPtrV[1] = 2 * colPtrM[1] + 1;
692 colPtrV[2] = 2 * colPtrM[2] + 1;
693
694 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 13), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
695
696 // FP = 21
697 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
698 valPtrM[0] = 0.140625;
699 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
700 valPtrM[1] = -0.046875;
701 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 2);
702 valPtrM[2] = 0.015625;
703 colPtrM[3] = (*coarseElementMDOFs)(coarseElement, 3);
704 valPtrM[3] = -0.046875;
705 colPtrM[4] = (*coarseElementMDOFs)(coarseElement, 4);
706 valPtrM[4] = 0.28125;
707 colPtrM[5] = (*coarseElementMDOFs)(coarseElement, 5);
708 valPtrM[5] = -0.09375;
709 colPtrM[6] = (*coarseElementMDOFs)(coarseElement, 6);
710 valPtrM[6] = -0.09375;
711 colPtrM[7] = (*coarseElementMDOFs)(coarseElement, 7);
712 valPtrM[7] = 0.28125;
713 colPtrM[8] = (*coarseElementMDOFs)(coarseElement, 8);
714 valPtrM[8] = 0.5625;
715
716 nnz = 9;
717 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[0], 8), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
718
719 colPtrV[0] = 2 * colPtrM[0];
720 colPtrV[1] = 2 * colPtrM[1];
721 colPtrV[2] = 2 * colPtrM[2];
722 colPtrV[3] = 2 * colPtrM[3];
723 colPtrV[4] = 2 * colPtrM[4];
724 colPtrV[5] = 2 * colPtrM[5];
725 colPtrV[6] = 2 * colPtrM[6];
726 colPtrV[7] = 2 * colPtrM[7];
727 colPtrV[8] = 2 * colPtrM[8];
728
729 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 16), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
730
731 colPtrV[0] = 2 * colPtrM[0] + 1;
732 colPtrV[1] = 2 * colPtrM[1] + 1;
733 colPtrV[2] = 2 * colPtrM[2] + 1;
734 colPtrV[3] = 2 * colPtrM[3] + 1;
735 colPtrV[4] = 2 * colPtrM[4] + 1;
736 colPtrV[5] = 2 * colPtrM[5] + 1;
737 colPtrV[6] = 2 * colPtrM[6] + 1;
738 colPtrV[7] = 2 * colPtrM[7] + 1;
739 colPtrV[8] = 2 * colPtrM[8] + 1;
740
741 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[0], 17), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
742
743 // FP = 22
744 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
745 valPtrM[0] = -0.046875;
746 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
747 valPtrM[1] = 0.140625;
748 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 2);
749 valPtrM[2] = -0.046875;
750 colPtrM[3] = (*coarseElementMDOFs)(coarseElement, 3);
751 valPtrM[3] = 0.015625;
752 colPtrM[4] = (*coarseElementMDOFs)(coarseElement, 4);
753 valPtrM[4] = 0.28125;
754 colPtrM[5] = (*coarseElementMDOFs)(coarseElement, 5);
755 valPtrM[5] = 0.28125;
756 colPtrM[6] = (*coarseElementMDOFs)(coarseElement, 6);
757 valPtrM[6] = -0.09375;
758 colPtrM[7] = (*coarseElementMDOFs)(coarseElement, 7);
759 valPtrM[7] = -0.09375;
760 colPtrM[8] = (*coarseElementMDOFs)(coarseElement, 8);
761 valPtrM[8] = 0.5625;
762
763 nnz = 9;
764 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[1], 8), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
765
766 colPtrV[0] = 2 * colPtrM[0];
767 colPtrV[1] = 2 * colPtrM[1];
768 colPtrV[2] = 2 * colPtrM[2];
769 colPtrV[3] = 2 * colPtrM[3];
770 colPtrV[4] = 2 * colPtrM[4];
771 colPtrV[5] = 2 * colPtrM[5];
772 colPtrV[6] = 2 * colPtrM[6];
773 colPtrV[7] = 2 * colPtrM[7];
774 colPtrV[8] = 2 * colPtrM[8];
775
776 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 16), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
777
778 colPtrV[0] = 2 * colPtrM[0] + 1;
779 colPtrV[1] = 2 * colPtrM[1] + 1;
780 colPtrV[2] = 2 * colPtrM[2] + 1;
781 colPtrV[3] = 2 * colPtrM[3] + 1;
782 colPtrV[4] = 2 * colPtrM[4] + 1;
783 colPtrV[5] = 2 * colPtrM[5] + 1;
784 colPtrV[6] = 2 * colPtrM[6] + 1;
785 colPtrV[7] = 2 * colPtrM[7] + 1;
786 colPtrV[8] = 2 * colPtrM[8] + 1;
787
788 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[1], 17), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
789
790 // FP = 23
791 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
792 valPtrM[0] = -0.046875;
793 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
794 valPtrM[1] = 0.015625;
795 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 2);
796 valPtrM[2] = -0.046875;
797 colPtrM[3] = (*coarseElementMDOFs)(coarseElement, 3);
798 valPtrM[3] = 0.140625;
799 colPtrM[4] = (*coarseElementMDOFs)(coarseElement, 4);
800 valPtrM[4] = -0.09375;
801 colPtrM[5] = (*coarseElementMDOFs)(coarseElement, 5);
802 valPtrM[5] = -0.09375;
803 colPtrM[6] = (*coarseElementMDOFs)(coarseElement, 6);
804 valPtrM[6] = 0.28125;
805 colPtrM[7] = (*coarseElementMDOFs)(coarseElement, 7);
806 valPtrM[7] = 0.28125;
807 colPtrM[8] = (*coarseElementMDOFs)(coarseElement, 8);
808 valPtrM[8] = 0.5625;
809
810 nnz = 9;
811 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[2], 8), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
812
813 colPtrV[0] = 2 * colPtrM[0];
814 colPtrV[1] = 2 * colPtrM[1];
815 colPtrV[2] = 2 * colPtrM[2];
816 colPtrV[3] = 2 * colPtrM[3];
817 colPtrV[4] = 2 * colPtrM[4];
818 colPtrV[5] = 2 * colPtrM[5];
819 colPtrV[6] = 2 * colPtrM[6];
820 colPtrV[7] = 2 * colPtrM[7];
821 colPtrV[8] = 2 * colPtrM[8];
822
823 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 16), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
824
825 colPtrV[0] = 2 * colPtrM[0] + 1;
826 colPtrV[1] = 2 * colPtrM[1] + 1;
827 colPtrV[2] = 2 * colPtrM[2] + 1;
828 colPtrV[3] = 2 * colPtrM[3] + 1;
829 colPtrV[4] = 2 * colPtrM[4] + 1;
830 colPtrV[5] = 2 * colPtrM[5] + 1;
831 colPtrV[6] = 2 * colPtrM[6] + 1;
832 colPtrV[7] = 2 * colPtrM[7] + 1;
833 colPtrV[8] = 2 * colPtrM[8] + 1;
834
835 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[2], 17), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
836
837 // FP = 24
838 colPtrM[0] = (*coarseElementMDOFs)(coarseElement, 0);
839 valPtrM[0] = 0.015625;
840 colPtrM[1] = (*coarseElementMDOFs)(coarseElement, 1);
841 valPtrM[1] = -0.046875;
842 colPtrM[2] = (*coarseElementMDOFs)(coarseElement, 2);
843 valPtrM[2] = 0.140625;
844 colPtrM[3] = (*coarseElementMDOFs)(coarseElement, 3);
845 valPtrM[3] = -0.046875;
846 colPtrM[4] = (*coarseElementMDOFs)(coarseElement, 4);
847 valPtrM[4] = -0.09375;
848 colPtrM[5] = (*coarseElementMDOFs)(coarseElement, 5);
849 valPtrM[5] = 0.28125;
850 colPtrM[6] = (*coarseElementMDOFs)(coarseElement, 6);
851 valPtrM[6] = 0.28125;
852 colPtrM[7] = (*coarseElementMDOFs)(coarseElement, 7);
853 valPtrM[7] = -0.09375;
854 colPtrM[8] = (*coarseElementMDOFs)(coarseElement, 8);
855 valPtrM[8] = 0.5625;
856
857 nnz = 9;
858 PM->insertGlobalValues((*fineElementMDOFs)(fineElement[3], 8), colPtrM.view(0, nnz), valPtrM.view(0, nnz));
859
860 colPtrV[0] = 2 * colPtrM[0];
861 colPtrV[1] = 2 * colPtrM[1];
862 colPtrV[2] = 2 * colPtrM[2];
863 colPtrV[3] = 2 * colPtrM[3];
864 colPtrV[4] = 2 * colPtrM[4];
865 colPtrV[5] = 2 * colPtrM[5];
866 colPtrV[6] = 2 * colPtrM[6];
867 colPtrV[7] = 2 * colPtrM[7];
868 colPtrV[8] = 2 * colPtrM[8];
869
870 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 16), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
871
872 colPtrV[0] = 2 * colPtrM[0] + 1;
873 colPtrV[1] = 2 * colPtrM[1] + 1;
874 colPtrV[2] = 2 * colPtrM[2] + 1;
875 colPtrV[3] = 2 * colPtrM[3] + 1;
876 colPtrV[4] = 2 * colPtrM[4] + 1;
877 colPtrV[5] = 2 * colPtrM[5] + 1;
878 colPtrV[6] = 2 * colPtrM[6] + 1;
879 colPtrV[7] = 2 * colPtrM[7] + 1;
880 colPtrV[8] = 2 * colPtrM[8] + 1;
881
882 PV->insertGlobalValues((*fineElementVDOFs)(fineElement[3], 17), colPtrV.view(0, nnz), valPtrM.view(0, nnz));
883
884 // FPr = 5
885 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 0);
886 valPtrP[0] = 0.25;
887 colPtrP[1] = (*coarseElementPDOFs)(coarseElement, 1);
888 valPtrP[1] = 0.25;
889 colPtrP[2] = (*coarseElementPDOFs)(coarseElement, 2);
890 valPtrP[2] = 0.25;
891 colPtrP[3] = (*coarseElementPDOFs)(coarseElement, 3);
892 valPtrP[3] = 0.25;
893
894 nnz = 4;
895 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[0], 2), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
896
897 // FPr = 6
898 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 1);
899 valPtrP[0] = 0.5;
900 colPtrP[1] = (*coarseElementPDOFs)(coarseElement, 2);
901 valPtrP[1] = 0.5;
902
903 nnz = 2;
904 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[1], 2), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
905
906 // FPr = 7
907 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 2);
908 valPtrP[0] = 1.0;
909
910 nnz = 1;
911 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[3], 2), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
912
913 // FPr = 8
914 colPtrP[0] = (*coarseElementPDOFs)(coarseElement, 2);
915 valPtrP[0] = 0.5;
916 colPtrP[1] = (*coarseElementPDOFs)(coarseElement, 3);
917 valPtrP[1] = 0.5;
918
919 nnz = 2;
920 PP->insertGlobalValues((*fineElementPDOFs)(fineElement[3], 3), colPtrP.view(0, nnz), valPtrP.view(0, nnz));
921
922 // Update counters:
923 if ((coarseElement + 1) % (nCoarseElements) == 0) // if the end of a row of c.g. elements
924 {
925 fineElement[0] = fineElement[3] + 1;
926 fineElement[1] = fineElement[0] + 1;
927 fineElement[2] = fineElement[0] + nFineElements;
928 fineElement[3] = fineElement[2] + 1;
929 } else {
930 fineElement[0] = fineElement[1] + 1;
931 fineElement[1] = fineElement[0] + 1;
932 fineElement[2] = fineElement[3] + 1;
933 fineElement[3] = fineElement[2] + 1;
934 }
935 } // END OF BUILD LOOP
936
937 // Loop over V rows
938 for (GO VRow = 0; VRow < fNV; VRow++) {
939 Teuchos::ArrayView<const LO> colPtr;
940 Teuchos::ArrayView<const SC> valPtr;
941
942 PV->getGlobalRowView(VRow, colPtr, valPtr);
943
944 // Can be directly inserted!
945 P->insertGlobalValues(VRow, colPtr, valPtr);
946 }
947
948 // Loop over P rows
949 for (GO PRow = 0; PRow < fNP; PRow++) {
950 Teuchos::ArrayView<const LO> colPtr;
951 Teuchos::ArrayView<const SC> valPtr;
952
953 // Now do pressure column:
954 PP->getGlobalRowView(PRow, colPtr, valPtr);
955
956 Teuchos::ArrayRCP<LO> newColPtr(colPtr.size(), nV);
957 for (LO jj = 0; jj < colPtr.size(); jj++) {
958 newColPtr[jj] += colPtr[jj];
959 }
960
961 // Insert into A
962 P->insertGlobalValues(PRow + fNV, newColPtr.view(0, colPtr.size()), valPtr);
963 }
964
965 // Loop over M rows
966 for (GO MRow = 0; MRow < fNM; MRow++) {
967 Teuchos::ArrayView<const LO> colPtr;
968 Teuchos::ArrayView<const SC> valPtr;
969
970 // Now do magnetics column:
971 PM->getGlobalRowView(MRow, colPtr, valPtr);
972
973 Teuchos::ArrayRCP<LO> newColPtr(colPtr.size(), nV + nP);
974 for (LO jj = 0; jj < colPtr.size(); jj++) {
975 newColPtr[jj] += colPtr[jj];
976 }
977
978 // Insert into A
979 P->insertGlobalValues(MRow + fNV + fNP, newColPtr.view(0, colPtr.size()), valPtr);
980 }
981
982 // Fill-complete all matrices
983 PV->fillComplete(colMapforPV, rowMapforPV);
984 PP->fillComplete(colMapforPP, rowMapforPP);
985 PM->fillComplete(colMapforPM, rowMapforPM);
986 P->fillComplete();
987
988 // Set prolongators on the coarse grid
989 Set(coarseLevel, "PV", PV);
990 Set(coarseLevel, "PP", PP);
991 Set(coarseLevel, "PM", PM);
992 Set(coarseLevel, "P", P);
993
994 Set(coarseLevel, "NV", nV);
995 Set(coarseLevel, "NP", nP);
996 Set(coarseLevel, "NM", nM);
997
998} // end buildp
999
1000} // namespace MueLu
1001
1002#define MUELU_GEOINTERPFACTORY_SHORT
1003#endif // MUELU_GEOINTERPFACTORY_DEF_HPP
MueLu_GeoInterpFactory_decl.hpp
MueLu_Level.hpp
MueLu::GeoInterpFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition MueLu_GeoInterpFactory_def.hpp:83
MueLu::GeoInterpFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Definition MueLu_GeoInterpFactory_def.hpp:42
MueLu::GeoInterpFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition MueLu_GeoInterpFactory_def.hpp:70
MueLu::Level
Class that holds all level-specific information.
Definition MueLu_Level.hpp:63
MueLu::Level::GetFactoryManager
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition MueLu_Level.cpp:73
MueLu::Level::Get
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
Definition MueLu_Level.hpp:154
MueLu
Namespace for MueLu classes and methods.
Definition MueLu_BrickAggregationFactory_decl.hpp:42
MueLu::Runtime1
@ Runtime1
Description of what is happening (more verbose)
Definition MueLu_VerbosityLevel.hpp:26