MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
29 namespace MueLu {
30 
31 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
32 GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GeoInterpFactory() {
33  GetOStream(Runtime1) << "I constructed a GeoInterpFactory object... Nothing else to do here." << std::endl;
34 }
35 
36 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
37 GeoInterpFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::~GeoInterpFactory() {
38  // Should be empty. All destruction should be handled by Level-based get stuff and RCP
39 }
40 
41 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
42 void 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 
69 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
70 void 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 
82 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
83 void 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::Get
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level-&gt;Get&lt; RCP&lt;Matrix&gt; &gt;(&quot;A&quot;, factory) if factory == NULL =&gt; use default factory.
Definition: MueLu_Level.hpp:154
Teuchos_SerialDenseMatrix.hpp
MueLu::GeoInterpFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeoInterpFactory_def.hpp:83
GO
GlobalOrdinal GO
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
Teuchos::ArrayView::size
size_type size() const
LO
LocalOrdinal LO
MueLu_Level.hpp
Xpetra_Matrix.hpp
Xpetra_MultiVector.hpp
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu::GeoInterpFactory::~GeoInterpFactory
virtual ~GeoInterpFactory()
Destructor.
Definition: MueLu_GeoInterpFactory_def.hpp:37
MueLu::Level::GetFactoryManager
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:73
Xpetra_VectorFactory.hpp
MueLu::Runtime1
Description of what is happening (more verbose)
Definition: MueLu_VerbosityLevel.hpp:26
Xpetra::MapFactory::createUniformContigMap
static Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > createUniformContigMap(UnderlyingLib lib, global_size_t numElements, const Teuchos::RCP< const Teuchos::Comm< int >> &comm)
Teuchos::ArrayRCP::view
ArrayView< T > view(size_type lowerOffset, size_type size) const