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