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