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