MueLu Version of the Day
Loading...
Searching...
No Matches
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
65#include <MueLu_Level.hpp>
66#include <MueLu_Monitor.hpp>
67#include <MueLu_Utilities.hpp>
70
71namespace MueLu {
72
73 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
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>
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>
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));
207 RCP<Matrix> PV = rcp(new CrsMatrixWrap(rowMapforPV,maxEntriesPerRowV));
208 RCP<Matrix> PP = rcp(new CrsMatrixWrap(rowMapforPP,maxEntriesPerRowP));
209 RCP<Matrix> PM = rcp(new CrsMatrixWrap(rowMapforPM,maxEntriesPerRowM));
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
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
virtual ~GeoInterpFactory()
Destructor.
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.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Definition: MueLu_Level.cpp:96
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
Namespace for MueLu classes and methods.
@ Runtime1
Description of what is happening (more verbose)