Finite Element Domain Decomposition Library
FEDDLib
Loading...
Searching...
No Matches
MeshStructured_def.hpp
1#ifndef MeshStructured_def_hpp
2#define MeshStructured_def_hpp
3
12
13namespace FEDD {
14template <class SC, class LO, class GO, class NO>
15MeshStructured<SC,LO,GO,NO>::MeshStructured():
16Mesh<SC,LO,GO,NO>(),
17coorRec(0),
18length(0),
19height(0),
20width(0)
21{ }
22
23template <class SC, class LO, class GO, class NO>
24MeshStructured<SC,LO,GO,NO>::MeshStructured(CommConstPtrConst_Type& comm):
25Mesh<SC,LO,GO,NO>(comm),
26coorRec(0),
27length(0),
28height(0),
29width(0)
30{ }
31
32template <class SC, class LO, class GO, class NO>
33MeshStructured<SC,LO,GO,NO>::~MeshStructured(){
34
35}
36
37template <class SC, class LO, class GO, class NO>
38void MeshStructured<SC,LO,GO,NO>::setGeometry2DRectangle(std::vector<double> coordinates, double l, double h){
39
40 coorRec = coordinates;
41 length = l;
42 height = h;
43
44 this->dim_ = 2;
45}
46
47template <class SC, class LO, class GO, class NO>
48void MeshStructured<SC,LO,GO,NO>::setGeometry3DBox(std::vector<double> coordinates, double l, double w, double h){
49
50 coorRec = coordinates;
51 length = l;
52 width = w;
53 height = h;
54
55 this->dim_ = 3;
56}
57
58template <class SC, class LO, class GO, class NO>
59void MeshStructured<SC,LO,GO,NO>::setRankRange(int numProcsCoarseSolve){
60 std::get<0>(this->rankRange_) = 0;
61 std::get<1>(this->rankRange_) = this->comm_->getSize() - 1 - numProcsCoarseSolve;
62}
63
64// template <class SC, class LO, class GO, class NO>
65// void MeshStructured<SC,LO,GO,NO>::buildMesh2DTPM(std::string FEType,
66// int N,
67// int M,
68// int numProcsCoarseSolve,
69// std::string underlyingLib){
70
71// buildMesh2D( FEType, N, M, numProcsCoarseSolve, underlyingLib );
72
73// setRankRange( numProcsCoarseSolve );
74
75// buildSurfaceLinesSquare();
76
77// }
78
79// template <class SC, class LO, class GO, class NO>
80// void MeshStructured<SC,LO,GO,NO>::buildMesh2DMiniTPM(std::string FEType,
81// int N,
82// int M,
83// int numProcsCoarseSolve,
84// std::string underlyingLib){
85
86// this->FEType_ = FEType;
87
88// this->numElementsGlob_ = 4;
89// int nmbPoints;
90
91// vec2D_int_ptr_Type elementsVec;
92// if (FEType=="P2") {
93// nmbPoints = 15;
94// elementsVec.reset(new std::vector<std::vector<int> >(this->numElementsGlob_,std::vector<int>(6,-1)));
95// } else if(FEType=="P1"){
96// nmbPoints = 6;
97// elementsVec.reset(new std::vector<std::vector<int> >(this->numElementsGlob_,std::vector<int>(3,-1)));
98// }
99
100// this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(2,0.0)));
101// this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
102// this->pointsUni_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(2,0.0)));
103// this->bcFlagUni_.reset(new std::vector<int> (nmbPoints,0));
104
105// Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
106// for (int i=0; i<nmbPoints; i++) {
107// pointsRepGlobMapping[i] = i;
108// }
109
110 // this->mapRepeated_.reset(new Map<LO,GO,NO>((GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
111
112// this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
113
114// double h = 0.1;
115// int counter = 0;
116// if (FEType=="P2") {
117// for (int i=0; i<3; i++) {
118// for (int j=0; j<5; j++) {
119// (*this->pointsRep_)[counter][0] = j*h/2;
120// (*this->pointsRep_)[counter][1] = i*h/2;
121
122// (*this->pointsUni_)[counter][0] = j*h/2;
123// (*this->pointsUni_)[counter][1] = i*h/2;
124// counter++;
125// }
126// }
127// } else if(FEType=="P1"){
128// for (int i=0; i<2; i++) {
129// for (int j=0; j<3; j++) {
130// (*this->pointsRep_)[counter][0] = j*h;
131// (*this->pointsRep_)[counter][1] = i*h;
132
133// (*this->pointsUni_)[counter][0] = j*h;
134// (*this->pointsUni_)[counter][1] = i*h;
135// counter++;
136// }
137// }
138// }
139
140// vec_int_ptr_Type elementFlag = Teuchos::rcp( new vec_int_Type( elementsVec->size(), 0 ) );
141
142// counter = 0;
143// int S=1;
144// int R=2;
145// int P2M = 2*(R+1)-1;
146// if (FEType=="P2") {
147
148// for (int s=0; s < S; s++) {
149// for (int r=0; r < R; r++) {
150
151// (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) ;
152// (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
153// (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
154
155// (*elementsVec)[counter][3] = 2*(r) +1 + 2*P2M * (s) ;
156// (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
157// (*elementsVec)[counter][5] = 2*(r+1) + 2*P2M * (s) +P2M ;
158
159// counter++;
160
161
162
163// (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s+1) ;
164// (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
165// (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
166
167// (*elementsVec)[counter][3] = 2*(r) + 2*P2M * (s) +P2M ;
168// (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
169// (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s+1) ;
170
171// counter++;
172// }
173// }
174
175// } else if(FEType=="P1") {
176
177// for (int s=0; s < S; s++) {
178// for (int r=0; r < R; r++) {
179
180// (*elementsVec)[counter][0] = r+1 + (R+1)* s;
181// (*elementsVec)[counter][1] = r + (R+1)* s;
182// (*elementsVec)[counter][2] = r+1 + (R+1) * (s+1);
183
184// counter++;
185
186// (*elementsVec)[counter][0] = r + (R+1) * (s+1);
187// (*elementsVec)[counter][1] = r + (R+1) * (s);
188// (*elementsVec)[counter][2] = r+1 + (R+1) * (s+1);
189
190// counter++;
191// }
192
193// }
194// }
195
196// setRankRange( numProcsCoarseSolve );
197
198// buildElementsClass( elementsVec, elementFlag );
199
200// buildSurfaceLinesSquareMiniTPM( FEType );
201
202// }
203
204
205template <class SC, class LO, class GO, class NO>
206void MeshStructured<SC,LO,GO,NO>::buildElementsClass( vec2D_int_ptr_Type elements, vec_int_ptr_Type elementFlag ){
207
208 this->elementsC_.reset(new Elements ( this->FEType_, this->dim_ ) );
209 bool setFlags = !elementFlag.is_null();
210 for (int i=0; i<elements->size(); i++) {
211 std::vector<LO> tmpElement;
212 for (int j=0; j<elements->at(i).size(); j++) {
213 tmpElement.push_back( (*elements)[i][j] );
214 }
215 if (setFlags){
216 FiniteElement fe( tmpElement, (*elementFlag)[i] );
217 this->elementsC_->addElement( fe );
218 }
219 else{
220 FiniteElement fe( tmpElement );
221 this->elementsC_->addElement( fe );
222 }
223 }
224}
225
226template <class SC, class LO, class GO, class NO>
227void MeshStructured<SC,LO,GO,NO>::buildSurfaceLinesSquare(){
228
229// for (int i=0; i<this->elementsC_->numberElements(); i++) {
230//
231// }
232
233}
234
235
236template <class SC, class LO, class GO, class NO>
237void MeshStructured<SC,LO,GO,NO>::buildMesh2D(std::string FEType,
238 int N,
239 int M,
240 int numProcsCoarseSolve){
241
242 using Teuchos::RCP;
243 using Teuchos::rcp;
244 using Teuchos::ScalarTraits;
245
246
247 TEUCHOS_TEST_FOR_EXCEPTION(!(M>=1),std::logic_error,"H/h is to small.");
248 TEUCHOS_TEST_FOR_EXCEPTION(this->comm_.is_null(),std::runtime_error,"comm_ is null.");
249
250 bool verbose (this->comm_->getRank() == 0);
251
252 setRankRange( numProcsCoarseSolve );
253
254 if (verbose) {
255 std::cout << std::endl;
256 }
257
258 int rank = this->comm_->getRank();
259 int size = this->comm_->getSize() - numProcsCoarseSolve;
260
261 SC h = length/(M*N);//add variable length/width/heigth
262 SC H = length/N;
263 GO nmbPoints_oneDir;
264
265 LO nmbElements;
266 LO nmbPoints;
267 vec2D_int_ptr_Type elementsVec;
268 vec_int_ptr_Type elementFlag;
269
270 if (FEType == "P0") {
271 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"implement P0.");
272 }
273 else if (FEType == "P1") {
274 nmbPoints_oneDir = N * (M+1) - (N-1);
275 nmbPoints = (M+1)*(M+1);
276 }
277 else if(FEType == "P2"){
278 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
279 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1);
280 }
281 else {
282 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong FE-Type, either P1 or P2.");
283 }
284
285 this->FEType_ = FEType;
286
287 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
288
289 this->numElementsGlob_ = 2*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
290
291 if (rank>=size) {
292 M = -1; // keine Schleife wird ausgefuehrt
293 nmbElements = 0;
294 nmbPoints = 0;
295 }
296 else{
297 nmbElements = 2*(M)*(M);
298 }
299
300 // P1 Mesh
301 if (FEType == "P1") {
302 if (verbose) {
303 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
304 }
305 if (verbose) {
306 std::cout << "-- Building P1 Points Repeated ... " << std::endl;
307 }
308
309 this->pointsRep_.reset(new vec2D_dbl_Type(nmbPoints,std::vector<double>(2,0.0)));
310 this->bcFlagRep_.reset(new vec_int_Type (nmbPoints,0));
311 elementsVec = Teuchos::rcp(new vec2D_int_Type(nmbElements,std::vector<int>(3,-1)));
312 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
313
314 int counter = 0;
315 int offset_x = (rank % N);
316 int offset_y = 0;
317
318 if ((rank % (N*N))>=N) {
319 offset_y = (int) (rank % (N*N))/(N);
320 }
321
322 for (int s=0; s < M+1; s++) {
323 for (int r=0; r < M+1; r++) {
324 (*this->pointsRep_)[counter][0] = r*h + offset_x * H;
325 if ((*this->pointsRep_)[counter][0]<100*ScalarTraits<SC>::eps() && (*this->pointsRep_)[counter][0]>-100*ScalarTraits<SC>::eps()) { (*this->pointsRep_)[counter][0]=0.0;}
326 (*this->pointsRep_)[counter][1] = s*h + offset_y * H;
327 if ((*this->pointsRep_)[counter][1]<100*ScalarTraits<SC>::eps() && (*this->pointsRep_)[counter][1]>-100*ScalarTraits<SC>::eps()) {(*this->pointsRep_)[counter][1]=0.0;}
328 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + offset_x*(M) + offset_y*(nmbPoints_oneDir)*M;
329 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-100*ScalarTraits<SC>::eps()) || (*this->pointsRep_)[counter][0] < (coorRec[0]+100*ScalarTraits<SC>::eps()) ||
330 (*this->pointsRep_)[counter][1] > (coorRec[1]+height-100*ScalarTraits<SC>::eps()) || (*this->pointsRep_)[counter][1] < (coorRec[1]+100*ScalarTraits<SC>::eps())) {
331
332 (*this->bcFlagRep_)[counter] = 1;
333 }
334 counter++;
335 }
336 }
337
338 if (verbose) {
339 std::cout << " done! --" << std::endl;
340 }
341
342 if (verbose) {
343 std::cout << "-- Building P1 Repeated and Unique Map ... " << std::flush;
344 }
345
346 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
347
348 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
349
350 if (verbose) {
351 std::cout << " done! --" << std::endl;
352 }
353
354 if (verbose) {
355 std::cout << "-- Building P1 Unique Points ... " << std::flush;
356 }
357
358 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(2,0.0)));
359 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
360 LO index;
361 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
362
363 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
364
365 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
366 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
367 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
368 }
369 if (verbose) {
370 std::cout << " done! --" << std::endl;
371 }
372
373
374 if (verbose) {
375 std::cout << "-- Building P1 Elements ... " << std::flush;
376 }
377 vec_int_ptr_Type elementFlag = Teuchos::rcp(new vec_int_Type( elementsVec->size(),0 ) );
378 counter = 0;
379 double x_ref, y_ref;
380
381 for (int s=0; s < M; s++) {
382 for (int r=0; r < M; r++) {
383
384 (*elementsVec)[counter][0] = r+1 + (M+1) * s;
385 (*elementsVec)[counter][1] = r + (M+1) * s ;
386 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
387 x_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(0) ) / 3.;
388 y_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(1) ) / 3.;
389 if ( x_ref>=0.3 && x_ref<=0.7) {
390 if ( y_ref>= 0.6) {
391 elementFlag->at(counter) = 1;
392 }
393 }
394
395 counter++;
396
397 (*elementsVec)[counter][0] = r + (M+1) * (s+1);
398 (*elementsVec)[counter][1] = r + (M+1) * (s);
399 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
400
401 x_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(0) ) / 3.;
402 y_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(1) ) / 3.;
403 if ( x_ref>=0.3 && x_ref<=0.7) {
404 if ( y_ref>= 0.6) {
405 elementFlag->at(counter) = 1;
406 }
407 }
408
409 counter++;
410 }
411 }
412
413 if (verbose) {
414 std::cout << " done! --" << std::endl;
415 }
416 }
417
418 // P2 Mesh
419
420
421 else if(FEType == "P2"){
422 if (verbose) {
423 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
424 }
425 if (verbose) {
426 std::cout << "-- Building P2 Points Repeated ... " << std::flush;
427 }
428
429 this->pointsRep_.reset(new vec2D_dbl_Type(nmbPoints, vec_dbl_Type(2,0.0)));
430 this->bcFlagRep_.reset(new vec_int_Type (nmbPoints,0));
431 elementsVec = Teuchos::rcp(new vec2D_int_Type(nmbElements, vec_int_Type(6,-1)));
432 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
433
434 int counter = 0;
435 int offset_x = (rank % N);
436 int offset_y = 0;
437
438 if ((rank % (N*N))>=N) {
439 offset_y = (int) (rank % (N*N))/(N);
440 }
441 bool p1point;
442 int p1_s = 0;
443 int p1_r = 0;
444 for (int s=0; s < 2*(M+1)-1; s++) {
445 for (int r=0; r < 2*(M+1)-1; r++) {
446 p1point = false;
447 if (s%2==0 && r%2==0) {
448 p1point = true;
449 p1_s = s/2;
450 p1_r = r/2;
451 }
452 (*this->pointsRep_)[counter][0] = r*h/2.0 + offset_x * H;
453 if ((*this->pointsRep_)[counter][0]<100*ScalarTraits<SC>::eps() && (*this->pointsRep_)[counter][0]>-100*ScalarTraits<SC>::eps()) (*this->pointsRep_)[counter][0]=0.0;
454 (*this->pointsRep_)[counter][1] = s*h/2.0 + offset_y * H;
455 if ((*this->pointsRep_)[counter][1]<100*ScalarTraits<SC>::eps() && (*this->pointsRep_)[counter][1]>-100*ScalarTraits<SC>::eps()) (*this->pointsRep_)[counter][1]=0.0;
456
457 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + offset_x*(2*(M+1)-2) + offset_y*(nmbPoints_oneDir)*(2*(M+1)-2) ;
458
459 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-100*ScalarTraits<SC>::eps()) || (*this->pointsRep_)[counter][0] < (coorRec[0]+100*ScalarTraits<SC>::eps()) ||
460 (*this->pointsRep_)[counter][1] > (coorRec[1]+height-100*ScalarTraits<SC>::eps()) || (*this->pointsRep_)[counter][1] < (coorRec[1]+100*ScalarTraits<SC>::eps())){
461 (*this->bcFlagRep_)[counter] = 1;
462
463 }
464 counter++;
465 }
466 }
467
468 if (verbose)
469 std::cout << " done! --" << std::endl;
470
471 if (verbose)
472 std::cout << "-- Building P2 Repeated and Unique Map ... " << std::flush;
473
474 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
475
476 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
477
478 if (verbose) {
479 std::cout << " done! --" << std::endl;
480 }
481
482 if (verbose) {
483 std::cout << "-- Building P2 Unique Points ... " << std::flush;
484 }
485
486 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(2,0.0)));
487 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
488
489 LO index;
490 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
491
492 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
493
494 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
495 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
496 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
497 }
498
499 if (verbose) {
500 std::cout << " done! --" << std::endl;
501 }
502
503 // Triangle numbering
504 // 2
505 // * *
506 // * *
507 // 4 5
508 // * *
509 // * *
510 // 1 * * 3 * * 0
511
512
513 if (verbose)
514 std::cout << "-- Building P2 Elements ... " << std::flush;
515
516 int P2M = 2*(M+1)-1;
517
518 vec_int_ptr_Type elementFlag = Teuchos::rcp(new vec_int_Type( elementsVec->size(),0 ) );
519 counter = 0;
520 double x_ref, y_ref;
521
522 for (int s=0; s < M; s++) {
523 for (int r=0; r < M; r++) {
524
525 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) ;
526 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
527 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
528
529 (*elementsVec)[counter][3] = 2*(r) +1 + 2*P2M * (s) ;
530 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
531 (*elementsVec)[counter][5] = 2*(r+1) + 2*P2M * (s) +P2M ;
532
533 x_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(0) ) / 3.;
534 y_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(1) ) / 3.;
535 if ( x_ref>=0.3 && x_ref<=0.7) {
536 if ( y_ref>= 0.6) {
537 elementFlag->at(counter) = 1;
538 }
539 }
540
541 counter++;
542
543 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s+1) ;
544 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
545 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
546
547 (*elementsVec)[counter][3] = 2*(r) + 2*P2M * (s) +P2M ;
548 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
549 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s+1) ;
550
551 x_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(0) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(0) ) / 3.;
552 y_ref = ( this->pointsRep_->at( elementsVec->at(counter).at(0) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(1) ).at(1) + this->pointsRep_->at( elementsVec->at(counter).at(2) ).at(1) ) / 3.;
553 if ( x_ref>=0.3 && x_ref<=0.7) {
554 if ( y_ref>= 0.6) {
555 elementFlag->at(counter) = 1;
556 }
557 }
558
559 counter++;
560
561 }
562 }
563
564
565
566 if (verbose) {
567 std::cout << " done! --" << std::endl;
568 }
569 }
570 buildElementsClass(elementsVec, elementFlag);
571
572}
573
574template <class SC, class LO, class GO, class NO>
575void MeshStructured<SC,LO,GO,NO>::buildMesh3D(std::string FEType,
576 int N,
577 int M,
578 int numProcsCoarseSolve){
579
580 using Teuchos::RCP;
581 using Teuchos::rcp;
582 using Teuchos::ScalarTraits;
583
584 TEUCHOS_TEST_FOR_EXCEPTION(!(M>=1),std::logic_error,"H/h is to small.");
585 TEUCHOS_TEST_FOR_EXCEPTION(this->comm_.is_null(),std::runtime_error,"comm_ is null.");
586
587 bool verbose (this->comm_->getRank() == 0);
588
589 setRankRange( numProcsCoarseSolve );
590
591 if (verbose) {
592 std::cout << std::endl;
593 }
594
595 SC eps = ScalarTraits<SC>::eps();
596
597 int rank = this->comm_->getRank();
598 int size = this->comm_->getSize() - numProcsCoarseSolve;
599
600 SC h = length/(M*N);//add variable length/width/heigth
601 SC H = length/N;
602
603 LO nmbPoints_oneDir;
604
605 LO nmbElements;
606 LO nmbPoints;
607 if (FEType == "P0") {
608 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"implement P0.");
609 }
610 else if (FEType == "P1") {
611 nmbPoints_oneDir = N * (M+1) - (N-1);
612 nmbPoints = (M+1)*(M+1)*(M+1);
613 }
614 else if (FEType == "P2"){
615 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
616 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
617 }
618 else if (FEType == "P2-CR"){
619 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"P2-CR might not work properly.");
620 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
621 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
622 }
623 else if (FEType == "P1-disc" || FEType == "P1-disc-global"){
624
625 }
626 else if (FEType == "Q1"){
627
628 }
629 else if (FEType == "Q2"){
630
631 }
632 else if (FEType == "Q2-20"){
633
634 }
635 else
636 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong FE-Type, only P1,P1-disc, P1-disc-global, P2, P2-CR, Q1, Q2, Q2-20.");
637
638 this->FEType_ = FEType;
639
640 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
641
642 this->numElementsGlob_ = 6*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
643 int MM=M;
644 if (rank>=size) {
645 M = -1; // keine Schleife wird ausgefuehrt
646 nmbElements = 0;
647 nmbPoints = 0;
648 }
649 else{
650 nmbElements = 6*M*M*M;
651 }
652 vec2D_int_ptr_Type elementsVec;
653 vec_int_ptr_Type elementFlag;
654 // P1 Mesh
655 if (FEType == "P1") {
656
657 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
658 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
659 elementsVec = Teuchos::rcp(new vec2D_int_Type( nmbElements, vec_int_Type(4, -1) ));
660 elementFlag = Teuchos::rcp(new vec_int_Type( elementsVec->size(),0 ) );
661 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
662
663 int counter = 0;
664 int offset_x = (rank % N);
665 int offset_y = 0;
666 int offset_z = 0;
667
668 if ((rank % (N*N))>=N) {
669 offset_y = (int) (rank % (N*N))/(N);
670 }
671
672 if ((rank % (N*N*N))>=N*N ) {
673 offset_z = (int) (rank % (N*N*N))/(N*(N));
674 }
675
676 for (int t=0; t < M+1; t++) {
677 for (int s=0; s < M+1; s++) {
678 for (int r=0; r < M+1; r++) {
679 (*this->pointsRep_)[counter][0] = r*h + offset_x * H;
680 if ((*this->pointsRep_)[counter][0]<eps && (*this->pointsRep_)[counter][0]>-eps) (*this->pointsRep_)[counter][0]=0.0;
681
682 (*this->pointsRep_)[counter][1] = s*h + offset_y * H;
683 if ((*this->pointsRep_)[counter][1]<eps && (*this->pointsRep_)[counter][1]>-eps) (*this->pointsRep_)[counter][1]=0.0;
684
685 (*this->pointsRep_)[counter][2] = t*h + offset_z * H;
686 if ((*this->pointsRep_)[counter][2]<eps && (*this->pointsRep_)[counter][2]>-eps) (*this->pointsRep_)[counter][2]=0.0;
687
688 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + t*nmbPoints_oneDir*nmbPoints_oneDir \
689 + offset_x*(M) + offset_y*(nmbPoints_oneDir)*M + offset_z*(nmbPoints_oneDir)*(nmbPoints_oneDir)*M ;
690
691 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
692 (*this->pointsRep_)[counter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps) ||
693 (*this->pointsRep_)[counter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[counter][2] < (coorRec[2]+eps) ) {
694
695 (*this->bcFlagRep_)[counter] = 1;
696 }
697
698 counter++;
699 }
700 }
701 }
702
703 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
704
705 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
706
707 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
708 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
709 LO index;
710 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
711
712 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
713
714 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
715 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
716 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
717 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
718 }
719
720 counter = 0;
721 for (int t=0; t < M; t++) {
722 for (int s=0; s < M; s++) {
723 for (int r=0; r < M; r++) {
724 (*elementsVec)[counter][0] = r+1 + (M+1) * s + (M+1)*(M+1) * t ;
725 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
726 (*elementsVec)[counter][2] = r+1 + (M+1) * s + (M+1)*(M+1) * (t+1) ;
727 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
728 counter++;
729 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * (t+1) ;
730 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
731 (*elementsVec)[counter][2] = r+1 + (M+1) * s + (M+1)*(M+1) * (t+1) ;
732 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
733 counter++;
734 (*elementsVec)[counter][0] = r+1 + (M+1) * s + (M+1)*(M+1) * t ;
735 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
736 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
737 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
738 counter++;
739 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
740 (*elementsVec)[counter][1] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
741 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
742 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
743 counter++;
744 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
745 (*elementsVec)[counter][1] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
746 (*elementsVec)[counter][2] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
747 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
748 counter++;
749 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
750 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * (t+1) ;
751 (*elementsVec)[counter][2] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
752 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
753 counter++;
754 }
755 }
756 }
757 buildElementsClass(elementsVec, elementFlag);
758 }
759
760 else if(FEType == "P2"){
761
762 this->pointsRep_.reset(new vec2D_dbl_Type(nmbPoints, vec_dbl_Type(3, 0.0)));
763 this->bcFlagRep_.reset(new vec_int_Type (nmbPoints, 0));
764 elementsVec = Teuchos::rcp(new vec2D_int_Type(nmbElements, vec_int_Type(10, -1)));
765 elementFlag = Teuchos::rcp(new vec_int_Type( elementsVec->size(),0 ) );
766 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
767
768 int counter = 0;
769 int offset_x = (rank % N);
770 int offset_y = 0;
771 int offset_z = 0;
772
773 if ((rank % (N*N))>=N) {
774 offset_y = (int) (rank % (N*N))/(N);
775 }
776
777 if ((rank % (N*N*N))>=N*N ) {
778 offset_z = (int) (rank % (N*N*N))/(N*(N));
779 }
780 bool p1point;
781 int p1_s = 0;
782 int p1_r = 0;
783 int p1_t = 0;
784 for (int t=0; t < 2*(M+1)-1; t++) {
785 for (int s=0; s < 2*(M+1)-1; s++) {
786 for (int r=0; r < 2*(M+1)-1; r++) {
787 p1point = false;
788 if (s%2==0 && r%2==0 && t%2==0) {
789 p1point = true;
790 p1_s = s/2;
791 p1_r = r/2;
792 p1_t = t/2;
793 }
794 (*this->pointsRep_)[counter][0] = r*h/2.0 + offset_x * H;
795 if ((*this->pointsRep_)[counter][0]<eps && (*this->pointsRep_)[counter][0]>-eps) (*this->pointsRep_)[counter][0]=0.0;
796 (*this->pointsRep_)[counter][1] = s*h/2.0 + offset_y * H;
797 if ((*this->pointsRep_)[counter][1]<eps && (*this->pointsRep_)[counter][1]>-eps) (*this->pointsRep_)[counter][1]=0.0;
798 (*this->pointsRep_)[counter][2] = t*h/2.0 + offset_z * H;
799 if ((*this->pointsRep_)[counter][2]<eps && (*this->pointsRep_)[counter][2]>-eps) (*this->pointsRep_)[counter][2]=0.0;
800
801 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + t*nmbPoints_oneDir*nmbPoints_oneDir \
802 + offset_x*(2*(M+1)-2) + offset_y*(nmbPoints_oneDir)*(2*(M+1)-2) + offset_z*(nmbPoints_oneDir)*(nmbPoints_oneDir)*(2*(M+1)-2) ;
803
804 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
805 (*this->pointsRep_)[counter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps) ||
806 (*this->pointsRep_)[counter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[counter][2] < (coorRec[2]+eps) ) {
807 (*this->bcFlagRep_)[counter] = 1;
808
809 }
810
811 counter++;
812 }
813 }
814 }
815
816 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
817
818 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
819
820 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
821 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
822
823 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
824 (*this->pointsUni_)[i][0] = (this->mapUnique_->getGlobalElement(i) % nmbPoints_oneDir) * h/2;
825 if ((*this->pointsUni_)[i][0]<eps && (*this->pointsUni_)[i][0]>-eps) (*this->pointsUni_)[i][0]=0.0;
826
827 (*this->pointsUni_)[i][1] = ((int) ((this->mapUnique_->getGlobalElement(i) % (nmbPoints_oneDir*nmbPoints_oneDir)) / nmbPoints_oneDir) + eps) *h/2;
828 if ((*this->pointsUni_)[i][1]<eps && (*this->pointsUni_)[i][1]>-eps) (*this->pointsUni_)[i][1]=0.0;
829
830 (*this->pointsUni_)[i][2] = ((int)(this->mapUnique_->getGlobalElement(i) / (nmbPoints_oneDir*nmbPoints_oneDir) + eps)) * h/2;
831 if ((*this->pointsUni_)[i][2]<eps && (*this->pointsUni_)[i][2]>-eps) (*this->pointsUni_)[i][2]=0.0;
832
833 if ((*this->pointsUni_)[i][0] > (coorRec[0]+length-eps) || (*this->pointsUni_)[i][0] < (coorRec[0]+eps) ||
834 (*this->pointsUni_)[i][1] > (coorRec[1]+width-eps) || (*this->pointsUni_)[i][1] < (coorRec[1]+eps) ||
835 (*this->pointsUni_)[i][2] > (coorRec[2]+height-eps) || (*this->pointsUni_)[i][2] < (coorRec[2]+eps) ) {
836 (*this->bcFlagUni_)[i] = 1;
837
838 }
839 }
840
841 // Face 1 Face2 Face 3 Face 4
842 // 2 2 * * 9 * * 3 3 * * 9 * * 2 3
843 // * * * * * * * *
844 // * * * * * * * *
845 // 5 6 6 7 8 5 8 7
846 // * * * * * * * *
847 // * * * * * * * *
848 // 1 * * 4 * * 0 0 1 1 * * 4 * * 0
849
850
851 int P2M = 2*(M+1)-1;
852
853 counter = 0;
854 for (int t=0; t < M; t++) {
855 for (int s=0; s < M; s++) {
856 for (int r=0; r < M; r++) {
857
858 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
859 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
860 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
861 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
862
863 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
864 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M;
865 (*elementsVec)[counter][7] = 2*(r+1) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M;
866 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M;
867 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M;
868 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1);
869
870 counter++;
871
872 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
873 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
874 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
875 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
876
877 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
878 (*elementsVec)[counter][6] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
879 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t+1) ;
880 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
881 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
882 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t+1);
883
884 counter++;
885
886 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
887 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
888 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
889 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
890
891 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
892 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) ;
893 (*elementsVec)[counter][7] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
894 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) ;
895 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
896 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
897
898 counter++;
899
900 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
901 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
902 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
903 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
904
905 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
906 (*elementsVec)[counter][6] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
907 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
908 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t);
909 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
910 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
911
912 counter++;
913
914 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
915 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
916 (*elementsVec)[counter][2] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
917 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
918
919 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
920 (*elementsVec)[counter][6] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
921 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
922 (*elementsVec)[counter][5] = 2*(r) + 2*P2M*(s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
923 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
924 (*elementsVec)[counter][9] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t+1) ;
925
926 counter++;
927
928 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
929 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
930 (*elementsVec)[counter][2] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
931 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
932
933 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
934 (*elementsVec)[counter][6] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
935 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
936 (*elementsVec)[counter][5] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) ;
937 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) ;
938 (*elementsVec)[counter][9] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t+1) ;
939
940 counter++;
941
942 }
943 }
944 }
945 buildElementsClass(elementsVec, elementFlag);
946 }
947 else if(FEType == "P1-disc" || FEType == "P1-disc-global")
948 buildP1_Disc_Q2_3DCube( N, MM, numProcsCoarseSolve );
949 else if(FEType == "Q1"){
950 build3DQ1Cube( N, M, numProcsCoarseSolve );
951 }
952 else if(FEType == "Q2"){
953 build3DQ2Cube( N, MM, numProcsCoarseSolve );
954 }
955 else if(FEType == "Q2-20"){
956 build3DQ2_20Cube( N, MM, numProcsCoarseSolve );
957 }
958
959
960
961};
962
963template <class SC, class LO, class GO, class NO>
964void MeshStructured<SC,LO,GO,NO>::buildP1_Disc_Q2_3DCube(int N,
965 int M,
966 int numProcsCoarseSolve){
967
968 using Teuchos::RCP;
969 using Teuchos::rcp;
970 using Teuchos::ScalarTraits;
971
972
973 bool verbose (this->comm_->getRank() == 0);
974
975 setRankRange( numProcsCoarseSolve );
976
977 if (verbose)
978 std::cout << std::endl;
979
980 SC eps = ScalarTraits<SC>::eps();
981
982 int rank = this->comm_->getRank();
983 int size = this->comm_->getSize() - numProcsCoarseSolve;
984
985 SC h = length/(M*N);//add variable length/width/heigth
986 SC H = length/N;
987
988 LO nmbPoints_oneDir;
989
990 LO nmbElements;
991 LO nmbPoints = 4*M*M*M; // 4 points for each element
992// nmbPoints_oneDir = N * (M+1) - (N-1);
993
994 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
995
996 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
997
998 if (rank>=size) {
999 M = -1; // keine Schleife wird ausgefuehrt
1000 nmbElements = 0;
1001 nmbPoints = 0;
1002 }
1003 else{
1004 nmbElements = M*M*M;
1005 }
1006
1007
1008 int counter = 0;
1009 int offset_x = (rank % N);
1010 int offset_y = 0;
1011 int offset_z = 0;
1012
1013 if ((rank % (N*N))>=N) {
1014 offset_y = (int) (rank % (N*N))/(N);
1015 }
1016
1017 if ((rank % (N*N*N))>=N*N ) {
1018 offset_z = (int) (rank % (N*N*N))/(N*(N));
1019 }
1020
1021
1022 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1023 this->pointsUni_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1024 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
1025 this->bcFlagUni_.reset(new std::vector<int> (nmbPoints,0));
1026 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
1027
1028
1029 if (verbose)
1030 std::cout << "-- Building P1-disc Points and Elements according to Q2 ... " << std::flush;
1031
1032 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new vec2D_int_Type(nmbElements, vec_int_Type(4, -1)));
1033
1034 counter = 0;
1035 LO pCounter = 0;
1036 GO globalCounterPoints = rank * nmbPoints;
1037 for (int t=0; t < M; t++) {
1038 for (int s=0; s < M; s++) {
1039 for (int r=0; r < M; r++) {
1040
1041 //point 1
1042 (*this->pointsRep_)[pCounter][0] = r*h + offset_x * H;
1043 (*this->pointsRep_)[pCounter][1] = s*h + offset_y * H;
1044 (*this->pointsRep_)[pCounter][2] = t*h + offset_z * H;
1045 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
1046 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
1047 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
1048
1049 (*this->bcFlagRep_)[counter] = 1;
1050 }
1051
1052 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
1053 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
1054 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
1055
1056 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
1057
1058 (*elementsVec)[counter][0] = pCounter;
1059 pointsRepGlobMapping[pCounter] = globalCounterPoints;
1060 globalCounterPoints++;
1061 pCounter++;
1062
1063 //point 2
1064 (*this->pointsRep_)[pCounter][0] = (r+1)*h + offset_x * H;
1065 (*this->pointsRep_)[pCounter][1] = s*h + offset_y * H;
1066 (*this->pointsRep_)[pCounter][2] = t*h + offset_z * H;
1067 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
1068 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
1069 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
1070
1071 (*this->bcFlagRep_)[counter] = 1;
1072 }
1073
1074 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
1075 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
1076 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
1077
1078 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
1079
1080 (*elementsVec)[counter][1] = pCounter;
1081 pointsRepGlobMapping[pCounter] = globalCounterPoints;
1082 globalCounterPoints++;
1083 pCounter++;
1084 //point 3
1085 (*this->pointsRep_)[pCounter][0] = r*h + offset_x * H;
1086 (*this->pointsRep_)[pCounter][1] = (s+1)*h + offset_y * H;
1087 (*this->pointsRep_)[pCounter][2] = t*h + offset_z * H;
1088 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
1089 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
1090 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
1091
1092 (*this->bcFlagRep_)[counter] = 1;
1093 }
1094
1095 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
1096 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
1097 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
1098
1099 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
1100
1101 (*elementsVec)[counter][2] = pCounter;
1102 pointsRepGlobMapping[pCounter] = globalCounterPoints;
1103 globalCounterPoints++;
1104 pCounter++;
1105
1106 //point 4
1107 (*this->pointsRep_)[pCounter][0] = r*h + offset_x * H;
1108 (*this->pointsRep_)[pCounter][1] = s*h + offset_y * H;
1109 (*this->pointsRep_)[pCounter][2] = (t+1)*h + offset_z * H;
1110 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
1111 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
1112 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
1113
1114 (*this->bcFlagRep_)[counter] = 1;
1115 }
1116
1117 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
1118 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
1119 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
1120
1121 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
1122
1123 (*elementsVec)[counter][3] = pCounter;
1124 pointsRepGlobMapping[pCounter] = globalCounterPoints;
1125 globalCounterPoints++;
1126 pCounter++;
1127
1128 counter++;
1129 }
1130 }
1131 }
1132 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1133
1134 this->mapUnique_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1135
1136 buildElementsClass(elementsVec);
1137
1138 if (verbose)
1139 std::cout << "done!" << std::endl;
1140
1141}
1142
1143template <class SC, class LO, class GO, class NO>
1144void MeshStructured<SC,LO,GO,NO>::build3DQ1Cube(int N,
1145 int M,
1146 int numProcsCoarseSolve)
1147{
1148
1149 using Teuchos::RCP;
1150 using Teuchos::rcp;
1151 using Teuchos::ScalarTraits;
1152
1153 bool verbose (this->comm_->getRank() == 0);
1154
1155 if (verbose)
1156 std::cout << std::endl;
1157
1158 setRankRange( numProcsCoarseSolve );
1159
1160 SC eps = ScalarTraits<SC>::eps();
1161
1162 int rank = this->comm_->getRank();
1163 int size = this->comm_->getSize() - numProcsCoarseSolve;
1164
1165 SC h = length/(M*N);//add variable length/width/heigth
1166 SC H = length/N;
1167
1168 LO nmbElements;
1169
1170 LO nmbPoints_oneDir = N * (M+1) - (N-1);
1171 LO nmbPoints = (M+1)*(M+1)*(M+1);
1172
1173 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
1174
1175 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
1176
1177 if (rank>=size) {
1178 M = -1; // keine Schleife wird ausgefuehrt
1179 nmbElements = 0;
1180 nmbPoints = 0;
1181 }
1182 else{
1183 nmbElements = M*M*M;
1184 }
1185
1186 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1187 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
1188 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements, std::vector<int>(8, -1)));
1189 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
1190
1191 int counter = 0;
1192 int offset_x = (rank % N);
1193 int offset_y = 0;
1194 int offset_z = 0;
1195
1196 if ((rank % (N*N))>=N) {
1197 offset_y = (int) (rank % (N*N))/(N);
1198 }
1199
1200 if ((rank % (N*N*N))>=N*N ) {
1201 offset_z = (int) (rank % (N*N*N))/(N*(N));
1202 }
1203 for (int t=0; t < M+1; t++) {
1204 for (int s=0; s < M+1; s++) {
1205 for (int r=0; r < M+1; r++) {
1206 (*this->pointsRep_)[counter][0] = r*h + offset_x * H;
1207 (*this->pointsRep_)[counter][1] = s*h + offset_y * H;
1208 (*this->pointsRep_)[counter][2] = t*h + offset_z * H;
1209
1210 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + t*nmbPoints_oneDir*nmbPoints_oneDir \
1211 + offset_x*(M) + offset_y*(nmbPoints_oneDir)*M + offset_z*(nmbPoints_oneDir)*(nmbPoints_oneDir)*M ;
1212
1213 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
1214 (*this->pointsRep_)[counter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps) ||
1215 (*this->pointsRep_)[counter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[counter][2] < (coorRec[2]+eps) ) {
1216
1217 (*this->bcFlagRep_)[counter] = 1;
1218 }
1219 counter++;
1220 }
1221 }
1222 }
1223 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1224
1225 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
1226
1227 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
1228 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
1229
1230 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
1231
1232 LO index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
1233
1234 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
1235 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
1236 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
1237 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
1238 }
1239
1240 LO offset = (M+1);
1241
1242 counter = 0;
1243 for (int t=0; t < M; t++) {
1244 for (int s=0; s < M; s++) {
1245 for (int r=0; r < M; r++) {
1246
1247 (*elementsVec)[counter][0] = r + (M+1) * (s) + (M+1)*(M+1) * t ;
1248 (*elementsVec)[counter][1] = r + 1 + (M+1) * (s) + (M+1)*(M+1) * t ;
1249 (*elementsVec)[counter][2] = r + 1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
1250 (*elementsVec)[counter][3] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
1251
1252 (*elementsVec)[counter][4] = r + (M+1) * (s) + (M+1)*(M+1) * (t+1) ;
1253 (*elementsVec)[counter][5] = r + 1 + (M+1) * (s) + (M+1)*(M+1) * (t+1) ;
1254 (*elementsVec)[counter][6] = r + 1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
1255 (*elementsVec)[counter][7] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
1256
1257 counter++;
1258
1259 }
1260 }
1261 }
1262 buildElementsClass(elementsVec);
1263}
1264
1265
1266template <class SC, class LO, class GO, class NO>
1267void MeshStructured<SC,LO,GO,NO>::build3DQ2Cube(int N,
1268 int M,
1269 int numProcsCoarseSolve)
1270{
1271
1272 using Teuchos::RCP;
1273 using Teuchos::rcp;
1274 using Teuchos::ScalarTraits;
1275
1276 bool verbose (this->comm_->getRank() == 0);
1277
1278 if (verbose)
1279 std::cout << std::endl;
1280
1281 setRankRange( numProcsCoarseSolve );
1282
1283 SC eps = ScalarTraits<SC>::eps();
1284
1285 int rank = this->comm_->getRank();
1286 int size = this->comm_->getSize() - numProcsCoarseSolve;
1287
1288 SC h = length/(M*N);//add variable length/width/heigth
1289 SC H = length/N;
1290
1291 LO nmbElements;
1292
1293 LO nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
1294 LO nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
1295
1296 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
1297
1298 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
1299
1300 if (rank>=size) {
1301 M = -1; // keine Schleife wird ausgefuehrt
1302 nmbElements = 0;
1303 nmbPoints = 0;
1304 }
1305 else{
1306 nmbElements = M*M*M;
1307 }
1308
1309 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1310 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
1311 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(27,-1)));
1312 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
1313
1314 int counter = 0;
1315 int offset_x = (rank % N);
1316 int offset_y = 0;
1317 int offset_z = 0;
1318
1319 if ((rank % (N*N))>=N) {
1320 offset_y = (int) (rank % (N*N))/(N);
1321 }
1322
1323 if ((rank % (N*N*N))>=N*N ) {
1324 offset_z = (int) (rank % (N*N*N))/(N*(N));
1325 }
1326 bool p1point;
1327 int p1_s = 0;
1328 int p1_r = 0;
1329 int p1_t = 0;
1330 for (int t=0; t < 2*(M+1)-1; t++) {
1331 for (int s=0; s < 2*(M+1)-1; s++) {
1332 for (int r=0; r < 2*(M+1)-1; r++) {
1333 p1point = false;
1334 if (s%2==0 && r%2==0 && t%2==0) {
1335 p1point = true;
1336 p1_s = s/2;
1337 p1_r = r/2;
1338 p1_t = t/2;
1339 }
1340 (*this->pointsRep_)[counter][0] = r*h/2.0 + offset_x * H;
1341 if ((*this->pointsRep_)[counter][0]<eps && (*this->pointsRep_)[counter][0]>-eps) (*this->pointsRep_)[counter][0]=0.0;
1342 (*this->pointsRep_)[counter][1] = s*h/2.0 + offset_y * H;
1343 if ((*this->pointsRep_)[counter][1]<eps && (*this->pointsRep_)[counter][1]>-eps) (*this->pointsRep_)[counter][1]=0.0;
1344 (*this->pointsRep_)[counter][2] = t*h/2.0 + offset_z * H;
1345 if ((*this->pointsRep_)[counter][2]<eps && (*this->pointsRep_)[counter][2]>-eps) (*this->pointsRep_)[counter][2]=0.0;
1346
1347 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + t*nmbPoints_oneDir*nmbPoints_oneDir \
1348 + offset_x*(2*(M+1)-2) + offset_y*(nmbPoints_oneDir)*(2*(M+1)-2) + offset_z*(nmbPoints_oneDir)*(nmbPoints_oneDir)*(2*(M+1)-2) ;
1349
1350 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
1351 (*this->pointsRep_)[counter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps) ||
1352 (*this->pointsRep_)[counter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[counter][2] < (coorRec[2]+eps) ) {
1353 (*this->bcFlagRep_)[counter] = 1;
1354
1355 }
1356
1357 counter++;
1358 }
1359 }
1360 }
1361
1362 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1363
1364 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
1365
1366 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
1367 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
1368
1369 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
1370 (*this->pointsUni_)[i][0] = (this->mapUnique_->getGlobalElement(i) % nmbPoints_oneDir) * h/2;
1371 if ((*this->pointsUni_)[i][0]<eps && (*this->pointsUni_)[i][0]>-eps) (*this->pointsUni_)[i][0]=0.0;
1372
1373 (*this->pointsUni_)[i][1] = ((int) ((this->mapUnique_->getGlobalElement(i) % (nmbPoints_oneDir*nmbPoints_oneDir)) / nmbPoints_oneDir) + eps) *h/2;
1374 if ((*this->pointsUni_)[i][1]<eps && (*this->pointsUni_)[i][1]>-eps) (*this->pointsUni_)[i][1]=0.0;
1375
1376 (*this->pointsUni_)[i][2] = ((int)(this->mapUnique_->getGlobalElement(i) / (nmbPoints_oneDir*nmbPoints_oneDir) + eps)) * h/2;
1377 if ((*this->pointsUni_)[i][2]<eps && (*this->pointsUni_)[i][2]>-eps) (*this->pointsUni_)[i][2]=0.0;
1378
1379 if ((*this->pointsUni_)[i][0] > (coorRec[0]+length-eps) || (*this->pointsUni_)[i][0] < (coorRec[0]+eps) ||
1380 (*this->pointsUni_)[i][1] > (coorRec[1]+width-eps) || (*this->pointsUni_)[i][1] < (coorRec[1]+eps) ||
1381 (*this->pointsUni_)[i][2] > (coorRec[2]+height-eps) || (*this->pointsUni_)[i][2] < (coorRec[2]+eps) ) {
1382 (*this->bcFlagUni_)[i] = 1;
1383
1384 }
1385 }
1386
1387 int P2M = 2*(M+1)-1;
1388
1389 counter = 0;
1390 for (int t=0; t < M; t++) {
1391 for (int s=0; s < M; s++) {
1392 for (int r=0; r < M; r++) {
1393
1394 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1395 (*elementsVec)[counter][1] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1396 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1397 (*elementsVec)[counter][3] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1398
1399 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1400 (*elementsVec)[counter][5] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1401 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1402 (*elementsVec)[counter][7] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1403
1404 (*elementsVec)[counter][8] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1405 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) ;
1406 (*elementsVec)[counter][10] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1407 (*elementsVec)[counter][11] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) ;
1408
1409 (*elementsVec)[counter][12] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1410 (*elementsVec)[counter][13] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1) ;
1411 (*elementsVec)[counter][14] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1412 (*elementsVec)[counter][15] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1) ;
1413
1414 (*elementsVec)[counter][16] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1415 (*elementsVec)[counter][17] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1416 (*elementsVec)[counter][18] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1417 (*elementsVec)[counter][19] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1418
1419 (*elementsVec)[counter][20] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1420 (*elementsVec)[counter][21] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1421 (*elementsVec)[counter][22] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1422 (*elementsVec)[counter][23] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1423
1424 (*elementsVec)[counter][24] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t);
1425 (*elementsVec)[counter][25] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1426 (*elementsVec)[counter][26] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1);
1427
1428 counter++;
1429 }
1430 }
1431 }
1432 buildElementsClass(elementsVec);
1433}
1434
1435template <class SC, class LO, class GO, class NO>
1436GO MeshStructured<SC,LO,GO,NO>::globalID_Q2_20Cube(int r, int s , int t, int &rr, int off_x, int off_y, int off_z, int M, int N, GO nmbPoints_oneDirFull, GO nmbPoints_oneDirMid){
1437
1438 GO index = -1;
1439 bool setPoint = false;
1440
1441 if (r%2==0 && t%2==0)
1442 setPoint = true;
1443 else{
1444 if (s%2==0 && t%2==0)
1445 setPoint = true;
1446 else{
1447 if (t%2==1 && r%2==0 && s%2==0) {
1448 setPoint = true;
1449 }
1450 }
1451 }
1452
1453
1454 if (setPoint) {
1455
1456 long long sizeFullSquare = nmbPoints_oneDirMid * nmbPoints_oneDirFull + nmbPoints_oneDirMid * (nmbPoints_oneDirMid-1);
1457 long long sizeNotFullSquare = nmbPoints_oneDirMid * nmbPoints_oneDirMid ;
1458 int ss = s/2;
1459 int tt = t/2;
1460
1461 index = rr + nmbPoints_oneDirFull * ( ss + (s%2) ) * (t%2==0) + nmbPoints_oneDirMid * ss;
1462 index += sizeFullSquare * ( tt + (t%2) ) + sizeNotFullSquare * tt;
1463 if (s%2==0)
1464 index += off_x * 2*M ;
1465 else
1466 index += off_x * M ;
1467
1468 index += off_y * ( nmbPoints_oneDirFull * M + nmbPoints_oneDirMid * M );
1469 index += off_z * ( sizeFullSquare * M + sizeNotFullSquare * M );
1470 rr++;
1471 }
1472
1473 return index;
1474}
1475
1476template <class SC, class LO, class GO, class NO>
1477void MeshStructured<SC,LO,GO,NO>::build3DQ2_20Cube(int N,
1478 int M,
1479 int numProcsCoarseSolve)
1480{
1481
1482 using Teuchos::RCP;
1483 using Teuchos::rcp;
1484 using Teuchos::ScalarTraits;
1485
1486 bool verbose (this->comm_->getRank() == 0);
1487
1488 setRankRange( numProcsCoarseSolve );
1489
1490 if (verbose)
1491 std::cout << std::endl;
1492 if (verbose)
1493 std::cout << "WARNING! Not working properly in parallel - fix global indexing." << std::endl;
1494
1495 SC eps = ScalarTraits<SC>::eps();
1496
1497 int rank = this->comm_->getRank();
1498 int size = this->comm_->getSize() - numProcsCoarseSolve;
1499
1500 SC h = length/(M*N);//add variable length/width/heigth
1501 SC H = length/N;
1502
1503 LO nmbElements;
1504
1505 LO nmbPoints_oneDirFull = N * (2*(M+1)-1) - (N-1);
1506 LO nmbPoints_oneDirMid = N * (M+1) - (N-1);
1507
1508 LO nmbPoints = ( (M+1) * (2*(M+1)-1) + M * (M+1) ) * (M+1) +
1509 ( (M+1) * (M+1) ) * M;
1510
1511 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
1512
1513 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
1514
1515 if (rank>=size) {
1516 M = -1; // keine Schleife wird ausgefuehrt
1517 nmbElements = 0;
1518 nmbPoints = 0;
1519 }
1520 else{
1521 nmbElements = M*M*M;
1522 }
1523
1524 this->pointsRep_.reset(new std::vector<std::vector<double> >(0));
1525 this->bcFlagRep_.reset(new std::vector<int> (0));
1526 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements, std::vector<int>(20,-1)));
1527 Teuchos::Array<GO> pointsRepGlobMapping(0);
1528
1529 int counter = 0;
1530 int offset_x = (rank % N);
1531 int offset_y = 0;
1532 int offset_z = 0;
1533
1534 if ((rank % (N*N))>=N) {
1535 offset_y = (int) (rank % (N*N))/(N);
1536 }
1537
1538 if ((rank % (N*N*N))>=N*N ) {
1539 offset_z = (int) (rank % (N*N*N))/(N*(N));
1540 }
1541 int rr=0;
1542 for (int t=0; t < 2*(M+1)-1; t++) {
1543 for (int s=0; s < 2*(M+1)-1; s++) {
1544 rr=0;
1545 for (int r=0; r < 2*(M+1)-1; r++) {
1546 GO index = globalID_Q2_20Cube( r, s, t, rr, offset_x, offset_y, offset_z, M, N,
1547 nmbPoints_oneDirFull, nmbPoints_oneDirMid );
1548
1549 if ( index>-1 ) {
1550 std::vector<double> p(3,0.0);
1551 p[0] = r*h/2.0 + offset_x * H;
1552 p[1] = s*h/2.0 + offset_y * H;
1553 p[2] = t*h/2.0 + offset_z * H;
1554 this->pointsRep_->push_back(p);
1555 pointsRepGlobMapping.push_back( index );
1556
1557 if (p[0] > (coorRec[0]+length-eps) || p[0] < (coorRec[0]+eps) ||
1558 p[1] > (coorRec[1]+width-eps) || p[1] < (coorRec[1]+eps) ||
1559 p[2] > (coorRec[2]+height-eps) || p[2] < (coorRec[2]+eps) )
1560 this->bcFlagRep_->push_back(1);
1561 else
1562 this->bcFlagRep_->push_back(0);
1563 counter++;
1564 }
1565 }
1566 }
1567 }
1568
1569 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1570
1571 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
1572
1573 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
1574 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
1575
1576 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
1577
1578 LO index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
1579
1580 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
1581 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
1582 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
1583 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
1584 }
1585
1586 int P2M = 2*(M+1)-1;
1587 int P1M = M+1;
1588 counter = 0;
1589 for (int t=0; t < M; t++) {
1590 for (int s=0; s < M; s++) {
1591 for (int r=0; r < M; r++) {
1592
1593 (*elementsVec)[counter][0] = 2*(r) + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1594 (*elementsVec)[counter][1] = 2*(r+1) + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1595 (*elementsVec)[counter][2] = 2*(r+1) + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1596 (*elementsVec)[counter][3] = 2*(r) + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1597
1598 (*elementsVec)[counter][4] = 2*(r) + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1599 (*elementsVec)[counter][5] = 2*(r+1) + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1600 (*elementsVec)[counter][6] = 2*(r+1) + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1601 (*elementsVec)[counter][7] = 2*(r) + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1602
1603 (*elementsVec)[counter][8] = 2*(r)+1 + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1604 (*elementsVec)[counter][9] = (r+1) + (P2M+P1M) * (s) + P2M + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1605 (*elementsVec)[counter][10] = 2*(r)+1 + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1606 (*elementsVec)[counter][11] = r + (P2M+P1M) * (s) + P2M + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) ;
1607
1608 (*elementsVec)[counter][12] = 2*(r)+1 + (P2M+P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1609 (*elementsVec)[counter][13] = (r+1) + (P2M+P1M) * (s) + P2M + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1610 (*elementsVec)[counter][14] = 2*(r)+1 + (P2M+P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1611 (*elementsVec)[counter][15] = r + (P2M+P1M) * (s) + P2M + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t+1) ;
1612
1613 (*elementsVec)[counter][16] = (r) + (P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) + P2M*(M+1)+P1M*M;
1614 (*elementsVec)[counter][17] = (r+1) + (P1M) * (s) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) + P2M*(M+1)+P1M*M;
1615 (*elementsVec)[counter][18] = (r+1) + (P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) + P2M*(M+1)+P1M*M;
1616 (*elementsVec)[counter][19] = r + (P1M) * (s+1) + ( (P2M*(M+1)+P1M*M) + (P1M*P1M) ) * (t) + P2M*(M+1)+P1M*M;
1617
1618 counter++;
1619 }
1620 }
1621 }
1622
1623}
1624
1625
1626template <class SC, class LO, class GO, class NO>
1627void MeshStructured<SC,LO,GO,NO>::build3DQ1BFS(int N,
1628 int M,
1629 int numProcsCoarseSolve){
1630
1631 using Teuchos::RCP;
1632 using Teuchos::rcp;
1633 using Teuchos::ScalarTraits;
1634
1635 typedef ScalarTraits<SC> ST;
1636 SC eps = ST::eps();
1637
1638 bool verbose (this->comm_->getRank() == 0);
1639
1640 setRankRange( numProcsCoarseSolve );
1641
1642 int rank = this->comm_->getRank();
1643 int size = this->comm_->getSize() - numProcsCoarseSolve;
1644
1645 int bfs_multiplier = (int) 2*(length)-1;
1646
1647 int nmbSubdomainsSquares = size / bfs_multiplier;
1648 int nmbSubdomainsSquares_OneDir = (std::pow(nmbSubdomainsSquares,1./3.) + 100*eps); // same as N
1649
1650 SC h = ST::one()/(M*N);
1651 SC H = ST::one()/N;
1652
1653 LO nmbPoints_oneDir = N * (M+1) - (N-1);
1654 LO nmbPoints = (M+1)*(M+1)*(M+1);
1655
1656
1657 //LO nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
1658 LO nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
1659 //LO nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
1660
1661
1662 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
1663 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1) * bfs_multiplier;
1664 LO nmbElements;
1665 if (rank>=size) {
1666 M = -1; // keine Schleife wird ausgefuehrt
1667 nmbElements = 0;
1668 nmbPoints = 0;
1669 }
1670 else{
1671 nmbElements = M*M*M;
1672 }
1673
1674 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1675 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,10));
1676 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements, std::vector<int>(27,-1)));
1677 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
1678
1679 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
1680
1681 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
1682 int offset_Squares_y = 0;
1683 int offset_Squares_z = ((whichSquareSet+1) % 2);
1684
1685 int counter = 0;
1686 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
1687 int offset_y = 0;
1688 int offset_z = 0;
1689 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
1690 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
1691 }
1692
1693 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))>=N*N ) {
1694 offset_z = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))/(N*(N));
1695 }
1696
1697 for (int t=0; t < M+1; t++) {
1698 for (int s=0; s < M+1; s++) {
1699 for (int r=0; r < M+1; r++) {
1700
1701 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
1702
1703 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
1704
1705 (*this->pointsRep_)[counter][2] = coorRec[2] + t*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
1706
1707 pointsRepGlobMapping[counter] = r
1708 + s*(nmbPoints_oneDir_allSubdomain);
1709 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1710 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
1711 }
1712 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
1713 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
1714 }
1715
1716 pointsRepGlobMapping[counter] += t*nmbPoints_oneDir_allSubdomain*nmbPoints_oneDir;
1717 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
1718 pointsRepGlobMapping[counter] -= t*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1719 }
1720
1721 pointsRepGlobMapping[counter] += offset_x*(M)
1722 + offset_y*( nmbPoints_oneDir_allSubdomain * M );
1723 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1724 pointsRepGlobMapping[counter] -= offset_y*M*(nmbPoints_oneDir-1);
1725 }
1726 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
1727 pointsRepGlobMapping[counter] -= offset_y*M*(nmbPoints_oneDir-1);
1728 }
1729
1730 pointsRepGlobMapping[counter] += offset_z * M * nmbPoints_oneDir_allSubdomain * nmbPoints_oneDir;
1731 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
1732 pointsRepGlobMapping[counter] -= offset_z*M*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1733 }
1734
1735 pointsRepGlobMapping[counter] += offset_Squares_x * M * nmbSubdomainsSquares_OneDir;
1736 if (offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1737 pointsRepGlobMapping[counter] -= M * nmbSubdomainsSquares_OneDir;
1738 }
1739 else if(offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
1740 pointsRepGlobMapping[counter] -= M * nmbSubdomainsSquares_OneDir;
1741 }
1742
1743 pointsRepGlobMapping[counter] += offset_Squares_z * nmbPoints_oneDir_allSubdomain * ((M) * nmbSubdomainsSquares_OneDir+1) * M * nmbSubdomainsSquares_OneDir;
1744 if (offset_Squares_z > 0 ) {
1745 pointsRepGlobMapping[counter] -= (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1746 }
1747 counter++;
1748 }
1749 }
1750 }
1751
1752 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1753
1754 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
1755
1756 if (verbose)
1757 std::cout << "-- Building Q2 Unique Points ... " << std::flush;
1758
1759 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
1760 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),10));
1761
1762 LO index;
1763 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
1764
1765 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
1766
1767 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
1768 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
1769 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
1770 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
1771 }
1772
1773 if (verbose)
1774 std::cout << " done! --" << std::endl;
1775
1776 //int P2M = 2*(M+1)-1;
1777
1778 counter = 0;
1779 for (int t=0; t < M; t++) {
1780 for (int s=0; s < M; s++) {
1781 for (int r=0; r < M; r++) {
1782
1783 (*elementsVec)[counter][0] = r + (M+1) * (s) + (M+1)*(M+1) * t ;
1784 (*elementsVec)[counter][1] = r + 1 + (M+1) * (s) + (M+1)*(M+1) * t ;
1785 (*elementsVec)[counter][2] = r + 1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
1786 (*elementsVec)[counter][3] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
1787
1788 (*elementsVec)[counter][4] = r + (M+1) * (s) + (M+1)*(M+1) * (t+1) ;
1789 (*elementsVec)[counter][5] = r + 1 + (M+1) * (s) + (M+1)*(M+1) * (t+1) ;
1790 (*elementsVec)[counter][6] = r + 1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
1791 (*elementsVec)[counter][7] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
1792
1793 counter++;
1794
1795 }
1796 }
1797 }
1798
1799 buildElementsClass(elementsVec);
1800}
1801
1802template <class SC, class LO, class GO, class NO>
1803void MeshStructured<SC,LO,GO,NO>::build3DQ2BFS(int N,
1804 int M,
1805 int numProcsCoarseSolve){
1806
1807 using Teuchos::RCP;
1808 using Teuchos::rcp;
1809 using Teuchos::ScalarTraits;
1810
1811 typedef ScalarTraits<SC> ST;
1812 SC eps = ST::eps();
1813
1814 bool verbose (this->comm_->getRank() == 0);
1815
1816 setRankRange( numProcsCoarseSolve );
1817
1818 int rank = this->comm_->getRank();
1819 int size = this->comm_->getSize() - numProcsCoarseSolve;
1820
1821 int bfs_multiplier = (int) 2*(length)-1;
1822
1823 int nmbSubdomainsSquares = size / bfs_multiplier;
1824 int nmbSubdomainsSquares_OneDir = (std::pow(nmbSubdomainsSquares,1./3.) + 100*eps); // same as N
1825
1826 SC h = ST::one()/(M*N);
1827 SC H = ST::one()/N;
1828
1829
1830 LO nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
1831 LO nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
1832 LO nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
1833
1834
1835 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
1836 this->numElementsGlob_ = (nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1) * bfs_multiplier;
1837 LO nmbElements;
1838 if (rank>=size) {
1839 M = -1; // keine Schleife wird ausgefuehrt
1840 nmbElements = 0;
1841 nmbPoints = 0;
1842 }
1843 else{
1844 nmbElements = M*M*M;
1845 }
1846
1847 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
1848 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
1849 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements, std::vector<int>(27,-1)));
1850 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
1851
1852 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
1853
1854 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
1855 int offset_Squares_y = 0;
1856 int offset_Squares_z = ((whichSquareSet+1) % 2);
1857
1858 int counter = 0;
1859 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
1860 int offset_y = 0;
1861 int offset_z = 0;
1862 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
1863 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
1864 }
1865
1866 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))>=N*N ) {
1867 offset_z = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))/(N*(N));
1868 }
1869
1870 for (int t=0; t < 2*(M+1)-1; t++) {
1871 for (int s=0; s < 2*(M+1)-1; s++) {
1872 for (int r=0; r < 2*(M+1)-1; r++) {
1873
1874 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h/2. + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
1875
1876 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h/2. + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
1877
1878 (*this->pointsRep_)[counter][2] = coorRec[2] + t*h/2. + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
1879
1880 pointsRepGlobMapping[counter] = r
1881 + s*(nmbPoints_oneDir_allSubdomain);
1882 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1883 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
1884 }
1885 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
1886 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
1887 }
1888
1889 pointsRepGlobMapping[counter] += t*nmbPoints_oneDir_allSubdomain*nmbPoints_oneDir;
1890 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
1891 pointsRepGlobMapping[counter] -= t*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1892 }
1893
1894 pointsRepGlobMapping[counter] += offset_x*(2*M)
1895 + offset_y*( nmbPoints_oneDir_allSubdomain * 2*M );
1896 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1897 pointsRepGlobMapping[counter] -= offset_y*2*M*(nmbPoints_oneDir-1);
1898 }
1899 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
1900 pointsRepGlobMapping[counter] -= offset_y*2*M*(nmbPoints_oneDir-1);
1901 }
1902
1903 pointsRepGlobMapping[counter] += offset_z * 2*M * nmbPoints_oneDir_allSubdomain * nmbPoints_oneDir;
1904 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
1905 pointsRepGlobMapping[counter] -= offset_z*2*M*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1906 }
1907
1908 pointsRepGlobMapping[counter] += offset_Squares_x * 2*M * nmbSubdomainsSquares_OneDir;
1909 if (offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
1910 pointsRepGlobMapping[counter] -= 2*M * nmbSubdomainsSquares_OneDir;
1911 }
1912 else if(offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
1913 pointsRepGlobMapping[counter] -= 2*M * nmbSubdomainsSquares_OneDir;
1914 }
1915
1916 pointsRepGlobMapping[counter] += offset_Squares_z * nmbPoints_oneDir_allSubdomain * ((2*M) * nmbSubdomainsSquares_OneDir+1) * 2*M * nmbSubdomainsSquares_OneDir;
1917 if (offset_Squares_z > 0 ) {
1918 pointsRepGlobMapping[counter] -= (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
1919 }
1920 counter++;
1921 }
1922 }
1923 }
1924
1925 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
1926
1927 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
1928
1929 if (verbose)
1930 std::cout << "-- Building Q2 Unique Points ... " << std::flush;
1931
1932 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
1933 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
1934
1935 LO index;
1936 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
1937
1938 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
1939
1940 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
1941 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
1942 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
1943 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
1944 }
1945
1946 if (verbose)
1947 std::cout << " done! --" << std::endl;
1948
1949 int P2M = 2*(M+1)-1;
1950
1951 counter = 0;
1952 for (int t=0; t < M; t++) {
1953 for (int s=0; s < M; s++) {
1954 for (int r=0; r < M; r++) {
1955
1956 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1957 (*elementsVec)[counter][1] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1958 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1959 (*elementsVec)[counter][3] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1960
1961 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1962 (*elementsVec)[counter][5] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1963 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1964 (*elementsVec)[counter][7] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1965
1966 (*elementsVec)[counter][8] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
1967 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) ;
1968 (*elementsVec)[counter][10] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
1969 (*elementsVec)[counter][11] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) ;
1970
1971 (*elementsVec)[counter][12] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
1972 (*elementsVec)[counter][13] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1) ;
1973 (*elementsVec)[counter][14] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
1974 (*elementsVec)[counter][15] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1) ;
1975
1976 (*elementsVec)[counter][16] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1977 (*elementsVec)[counter][17] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1978 (*elementsVec)[counter][18] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1979 (*elementsVec)[counter][19] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1980
1981 (*elementsVec)[counter][20] = 2*(r)+1 + 2*P2M * (s) + 2*P2M*P2M * (t) + P2M*P2M;
1982 (*elementsVec)[counter][21] = 2*(r+1) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1983 (*elementsVec)[counter][22] = 2*(r)+1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) + P2M*P2M;
1984 (*elementsVec)[counter][23] = 2*(r) + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1985
1986 (*elementsVec)[counter][24] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t);
1987 (*elementsVec)[counter][25] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t) + P2M*P2M;
1988 (*elementsVec)[counter][26] = 2*(r)+1 + 2*P2M * (s) + P2M + 2*P2M*P2M * (t+1);
1989
1990 counter++;
1991
1992 }
1993 }
1994 }
1995 buildElementsClass(elementsVec);
1996}
1997
1998template <class SC, class LO, class GO, class NO>
1999void MeshStructured<SC,LO,GO,NO>::buildMesh2DBFS(std::string FEType,
2000 int N,
2001 int M,
2002 int numProcsCoarseSolve) {
2003
2004
2005 using Teuchos::RCP;
2006 using Teuchos::rcp;
2007 using Teuchos::ScalarTraits;
2008
2009 TEUCHOS_TEST_FOR_EXCEPTION(!(M>=1),std::logic_error,"H/h is to small.");
2010 TEUCHOS_TEST_FOR_EXCEPTION(this->comm_.is_null(),std::runtime_error,"comm_ is null.");
2011
2012 SC eps = ScalarTraits<SC>::eps();
2013
2014 bool verbose (this->comm_->getRank() == 0);
2015
2016 setRankRange( numProcsCoarseSolve );
2017
2018 int rank = this->comm_->getRank();
2019 int size = this->comm_->getSize() - numProcsCoarseSolve;
2020
2021 int bfs_multiplier = (int) 2*(length)-1;
2022
2023 int nmbSubdomainsSquares = size / bfs_multiplier;
2024 int nmbSubdomainsSquares_OneDir = (std::pow(nmbSubdomainsSquares,1/2.) + 100*eps);
2025
2026 vec2D_int_ptr_Type elementsVec;
2027
2028 LO nmbElements;
2029 LO nmbPoints;
2030
2031 double h = 1./(M*N);
2032 double H = 1./N;
2033 GO nmbPoints_oneDir;
2034 GO nmbPoints_oneDir_allSubdomain;
2035 if (FEType == "P0") {
2036 nmbPoints_oneDir = N * (M+1) - (N-1) ;
2037 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
2038 nmbPoints = (M+1)*(M+1) ;
2039 }
2040 else if (FEType == "P1") {
2041 nmbPoints_oneDir = N * (M+1) - (N-1) ;
2042 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
2043 nmbPoints = (M+1)*(M+1) ;
2044 }
2045 else if(FEType == "P2"){
2046 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1) ;
2047 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
2048 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1) ;
2049 }
2050 else {
2051 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong FE-Type, either P1 or P2.");
2052 }
2053
2054 this->FEType_ = FEType;
2055
2056 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
2057
2058 this->numElementsGlob_ = 2*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1) * bfs_multiplier;
2059
2060 if (rank>=size) {
2061 M = -1; // keine Schleife wird ausgefuehrt
2062 nmbElements = 0;
2063 nmbPoints = 0;
2064 }
2065 else{
2066 nmbElements = 2*(M)*(M);
2067 }
2068
2069 // P0 Mesh
2070 if (FEType == "P0") {
2071 if (verbose)
2072 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
2073
2074 if (verbose)
2075 std::cout << "-- Building P0 Points Repeated ... " << std::endl;
2076
2077 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(2,0.0)));
2078 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2079 elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements, std::vector<int>(3,-1)));
2080
2081 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2082
2083 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2084
2085 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2086 int offset_Squares_y = ((whichSquareSet+1) % 2);
2087 int counter = 0;
2088 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2089 int offset_y = 0;
2090
2091 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
2092 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2093 }
2094
2095 for (int s=0; s < M+1; s++) {
2096 for (int r=0; r < M+1; r++) {
2097 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
2098 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
2099 pointsRepGlobMapping[counter] = r + s*(nmbPoints_oneDir_allSubdomain - (1-offset_Squares_y)*(nmbPoints_oneDir-1))
2100 + offset_x*(M)
2101 + offset_y*((nmbPoints_oneDir_allSubdomain) - (1-offset_Squares_y)*(nmbPoints_oneDir-1)) *M
2102 + offset_Squares_x * M * nmbSubdomainsSquares_OneDir
2103 + offset_Squares_y * (nmbPoints_oneDir_allSubdomain * M * nmbSubdomainsSquares_OneDir
2104 - (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1));
2105 //- M * nmbSubdomainsSquares_OneDir);
2106 if (offset_Squares_x>0 && offset_Squares_y==0 ) {
2107 pointsRepGlobMapping[counter] -= nmbPoints_oneDir-1;
2108 }
2109 if (offset_Squares_x>0 && offset_Squares_y==0 && offset_y+1==nmbSubdomainsSquares_OneDir && s==M) {
2110 pointsRepGlobMapping[counter] += nmbPoints_oneDir-1;
2111 }
2112 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
2113 (*this->pointsRep_)[counter][1] > (coorRec[1]+height-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps)) {
2114 (*this->bcFlagRep_)[counter] = 1;
2115 }
2116
2117 counter++;
2118 }
2119 }
2120
2121 if (verbose) {
2122 std::cout << " done! --" << std::endl;
2123 }
2124
2125 if (verbose) {
2126 std::cout << "-- Building P0 Repeated and Unique Map ... " << std::flush;
2127 }
2128
2129
2130
2131 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
2132
2133 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
2134 if (verbose) {
2135 std::cout << " done! --" << std::endl;
2136 }
2137
2138 if (verbose) {
2139 std::cout << "-- Building P0 Unique Points ... " << std::flush;
2140 }
2141
2142 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(), std::vector<double>(2,0.0)));
2143 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
2144 LO index;
2145 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
2146
2147 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
2148
2149 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
2150 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
2151 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
2152 }
2153
2154 if (verbose) {
2155 std::cout << " done! --" << std::endl;
2156 }
2157
2158
2159 if (verbose) {
2160 std::cout << "-- Building P0 Elements ... " << std::flush;
2161 }
2162
2163 counter = 0;
2164 for (int s=0; s < M; s++) {
2165 for (int r=0; r < M; r++) {
2166 (*elementsVec)[counter][0] = r+1 + (M+1) * s;
2167 (*elementsVec)[counter][1] = r + (M+1) * s ;
2168 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
2169 counter++;
2170 (*elementsVec)[counter][0] = r + (M+1) * (s+1);
2171 (*elementsVec)[counter][1] = r + (M+1) * (s);
2172 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
2173 counter++;
2174 }
2175 }
2176
2177 if (verbose) {
2178 std::cout << " done! --" << std::endl;
2179 }
2180 }
2181
2182 // P1 Mesh
2183 else if (FEType == "P1") {
2184 if (verbose) {
2185 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
2186 }
2187 if (verbose) {
2188 std::cout << "-- Building P1 Points Repeated ... " << std::endl;
2189 }
2190
2191 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(2,0.0)));
2192 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2193 elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(3,-1)));
2194
2195 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2196
2197 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2198
2199 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2200 int offset_Squares_y = ((whichSquareSet+1) % 2);
2201 int counter = 0;
2202 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2203 int offset_y = 0;
2204
2205 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
2206 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2207 }
2208
2209 for (int s=0; s < M+1; s++) {
2210 for (int r=0; r < M+1; r++) {
2211 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
2212 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
2213 pointsRepGlobMapping[counter] = r + s*(nmbPoints_oneDir_allSubdomain - (1-offset_Squares_y)*(nmbPoints_oneDir-1))
2214 + offset_x*(M)
2215 + offset_y*((nmbPoints_oneDir_allSubdomain) - (1-offset_Squares_y)*(nmbPoints_oneDir-1)) *M
2216 + offset_Squares_x * M * nmbSubdomainsSquares_OneDir
2217 + offset_Squares_y * (nmbPoints_oneDir_allSubdomain * M * nmbSubdomainsSquares_OneDir
2218 - (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1));
2219 //- M * nmbSubdomainsSquares_OneDir);
2220 if (offset_Squares_x>0 && offset_Squares_y==0 ) {
2221 pointsRepGlobMapping[counter] -= nmbPoints_oneDir-1;
2222 }
2223 if (offset_Squares_x>0 && offset_Squares_y==0 && offset_y+1==nmbSubdomainsSquares_OneDir && s==M) {
2224 pointsRepGlobMapping[counter] += nmbPoints_oneDir-1;
2225 }
2226 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
2227 (*this->pointsRep_)[counter][1] > (coorRec[1]+height-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps)) {
2228 (*this->bcFlagRep_)[counter] = 1;
2229 }
2230
2231 counter++;
2232 }
2233 }
2234
2235 if (verbose) {
2236 std::cout << " done! --" << std::endl;
2237 }
2238
2239 if (verbose) {
2240 std::cout << "-- Building P1 Repeated and Unique Map ... " << std::flush;
2241 }
2242
2243
2244
2245
2246 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
2247
2248 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
2249
2250 if (verbose) {
2251 std::cout << " done! --" << std::endl;
2252 }
2253
2254 if (verbose) {
2255 std::cout << "-- Building P1 Unique Points ... " << std::flush;
2256 }
2257
2258 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(), std::vector<double>(2,0.0)));
2259 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
2260 LO index;
2261 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
2262
2263 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
2264
2265 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
2266 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
2267 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
2268 }
2269
2270 if (verbose) {
2271 std::cout << " done! --" << std::endl;
2272 }
2273
2274
2275 if (verbose) {
2276 std::cout << "-- Building P1 Elements ... " << std::flush;
2277 }
2278
2279 counter = 0;
2280 for (int s=0; s < M; s++) {
2281 for (int r=0; r < M; r++) {
2282 (*elementsVec)[counter][0] = r+1 + (M+1) * s;
2283 (*elementsVec)[counter][1] = r + (M+1) * s ;
2284 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
2285 counter++;
2286 (*elementsVec)[counter][0] = r + (M+1) * (s+1);
2287 (*elementsVec)[counter][1] = r + (M+1) * (s);
2288 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1);
2289 counter++;
2290 }
2291 }
2292
2293 if (verbose) {
2294 std::cout << " done! --" << std::endl;
2295 }
2296 }
2297
2298 // P2 Mesh
2299 else if(FEType == "P2"){
2300 if (verbose) {
2301 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
2302 }
2303 if (verbose) {
2304 std::cout << "-- Building P2 Points Repeated ... " << std::flush;
2305 }
2306
2307 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(2,0.0)));
2308 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2309 elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(6,-1)));
2310 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2311
2312 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2313
2314 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2315 int offset_Squares_y = ((whichSquareSet+1) % 2);
2316
2317 int counter = 0;
2318 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2319 int offset_y = 0;
2320
2321
2322 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
2323 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2324 }
2325
2326
2327 bool p1point;
2328 int p1_s = 0;
2329 int p1_r = 0;
2330
2331 for (int s=0; s < 2*(M+1)-1; s++) {
2332 for (int r=0; r < 2*(M+1)-1; r++) {
2333 p1point = false;
2334 if (s%2==0 && r%2==0) {
2335 p1point = true;
2336 p1_s = s/2;
2337 p1_r = r/2;
2338 }
2339 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h/2. + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
2340 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h/2. + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
2341 pointsRepGlobMapping[counter] = r + s*(nmbPoints_oneDir_allSubdomain - (1-offset_Squares_y)*(nmbPoints_oneDir-1))
2342 + offset_x*(2*(M+1)-2)
2343 + offset_y*((nmbPoints_oneDir_allSubdomain) - (1-offset_Squares_y)*(nmbPoints_oneDir-1)) * (2*(M+1)-2)
2344 + offset_Squares_x * (2*(M+1)-2) * nmbSubdomainsSquares_OneDir
2345 + offset_Squares_y * (nmbPoints_oneDir_allSubdomain * 2*M * nmbSubdomainsSquares_OneDir
2346 - (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1));
2347 //- M * nmbSubdomainsSquares_OneDir);
2348 if (offset_Squares_x>0 && offset_Squares_y==0 ) {
2349 pointsRepGlobMapping[counter] -= nmbPoints_oneDir-1;
2350 }
2351 if (offset_Squares_x>0 && offset_Squares_y==0 && offset_y+1==nmbSubdomainsSquares_OneDir && s==2*M) {
2352 pointsRepGlobMapping[counter] += nmbPoints_oneDir-1;
2353 }
2354 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
2355 (*this->pointsRep_)[counter][1] > (coorRec[1]+height-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps)) {
2356 (*this->bcFlagRep_)[counter] = 1;
2357 }
2358
2359 counter++;
2360 }
2361 }
2362
2363 if (verbose) {
2364 std::cout << " done! --" << std::endl;
2365 }
2366 // int* globIndex = new int[MappingPointsRepGlob->size()];
2367 // globIndex = &(MappingPointsRepGlob->at(0));
2368
2369 if (verbose) {
2370 std::cout << "-- Building P2 Repeated and Unique Map ... " << std::flush;
2371 }
2372
2373
2374 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
2375
2376 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
2377
2378
2379 if (verbose) {
2380 std::cout << " done! --" << std::endl;
2381 }
2382
2383 if (verbose) {
2384 std::cout << "-- Building P2 Unique Points ... " << std::flush;
2385 }
2386
2387 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(2,0.0)));
2388 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
2389
2390 LO index;
2391 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
2392
2393 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
2394
2395 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
2396 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
2397 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
2398 }
2399
2400 if (verbose) {
2401 std::cout << " done! --" << std::endl;
2402 }
2403
2404 // Triangle numbering
2405 // 2
2406 // * *
2407 // * *
2408 // 4 5
2409 // * *
2410 // * *
2411 // 1 * * 3 * * 0
2412
2413
2414 if (verbose) {
2415 std::cout << "-- Building P2 Elements ... " << std::flush;
2416 }
2417
2418 int P2M = 2*(M+1)-1;
2419
2420 counter = 0;
2421
2422 for (int s=0; s < M; s++) {
2423 for (int r=0; r < M; r++) {
2424
2425 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) ;
2426 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
2427 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
2428
2429 (*elementsVec)[counter][3] = 2*(r) +1 + 2*P2M * (s) ;
2430 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
2431 (*elementsVec)[counter][5] = 2*(r+1) + 2*P2M * (s) +P2M ;
2432
2433 counter++;
2434
2435 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s+1) ;
2436 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) ;
2437 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) ;
2438
2439 (*elementsVec)[counter][3] = 2*(r) + 2*P2M * (s) +P2M ;
2440 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) +P2M ;
2441 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s+1) ;
2442
2443 counter++;
2444
2445 }
2446 }
2447
2448 if (verbose) {
2449 std::cout << " done! --" << std::endl;
2450 }
2451 }
2452 buildElementsClass(elementsVec);
2453}
2454
2455template <class SC, class LO, class GO, class NO>
2456void MeshStructured<SC,LO,GO,NO>::buildMesh3DBFS(std::string FEType,
2457 int N,
2458 int M,
2459 int numProcsCoarseSolve){
2460
2461 using Teuchos::RCP;
2462 using Teuchos::rcp;
2463 using Teuchos::ScalarTraits;
2464
2465 TEUCHOS_TEST_FOR_EXCEPTION(!(M>=1),std::logic_error,"H/h is to small.");
2466 TEUCHOS_TEST_FOR_EXCEPTION(this->comm_.is_null(),std::runtime_error,"comm_ is null.");
2467
2468 typedef ScalarTraits<SC> ST;
2469 SC eps = ST::eps();
2470
2471 bool verbose (this->comm_->getRank() == 0);
2472
2473 setRankRange( numProcsCoarseSolve );
2474
2475 int rank = this->comm_->getRank();
2476 int size = this->comm_->getSize() - numProcsCoarseSolve;
2477
2478 int bfs_multiplier = (int) 2*(length)-1;
2479
2480 int nmbSubdomainsSquares = size / bfs_multiplier;
2481 int nmbSubdomainsSquares_OneDir = (std::pow(nmbSubdomainsSquares,1./3.) + 100*eps); // same as N
2482
2483 SC h = ST::one()/(M*N);
2484 SC H = ST::one()/N;
2485
2486 LO nmbElements;
2487 LO nmbPoints;
2488
2489 GO nmbPoints_oneDir;
2490 GO nmbPoints_oneDir_allSubdomain;
2491
2492 if (FEType == "P0") {
2493 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"implement P0.");
2494 }
2495 else if (FEType == "P1") {
2496 nmbPoints_oneDir = N * (M+1) - (N-1);
2497 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1);
2498 nmbPoints = (M+1)*(M+1)*(M+1);
2499 }
2500 else if(FEType == "P2"){
2501 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
2502 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
2503 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
2504 }
2505 else if(FEType == "P2-CR"){
2506 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"P2-CR might not work properly.");
2507 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
2508 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1) ;
2509 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1);
2510 }
2511 else if(FEType == "P1-disc" || FEType == "P1-disc-global"){
2512 nmbPoints_oneDir = N * (M+1) - (N-1);
2513 nmbPoints_oneDir_allSubdomain = length * nmbPoints_oneDir - (length-1);
2514 nmbPoints = (M+1)*(M+1)*(M+1);
2515 }
2516 else if(FEType == "Q1"){
2517
2518 }
2519 else if(FEType == "Q2"){
2520
2521 }
2522 else
2523 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong FE-Type, only P1,P1-disc, P1-disc-global, P2, or P2-CR.");
2524
2525 this->FEType_ = FEType;
2526
2527 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
2528 this->numElementsGlob_ = 6*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1) * bfs_multiplier;
2529 int MM=M;
2530 if (rank>=size) {
2531 M = -1; // keine Schleife wird ausgefuehrt
2532 nmbElements = 0;
2533 nmbPoints = 0;
2534 }
2535 else{
2536 nmbElements = 6*M*M*M;
2537 }
2538
2539 // P1 Mesh
2540 if (FEType == "P1") {
2541
2542 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
2543 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2544 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(4,-1)));
2545
2546 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2547
2548 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2549
2550 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2551 int offset_Squares_y = 0;
2552 int offset_Squares_z = ((whichSquareSet+1) % 2);
2553
2554 int counter = 0;
2555 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2556 int offset_y = 0;
2557 int offset_z = 0;
2558 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N)
2559 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2560 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))>=N*N )
2561 offset_z = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))/(N*(N));
2562
2563 for (int t=0; t < M+1; t++) {
2564 for (int s=0; s < M+1; s++) {
2565 for (int r=0; r < M+1; r++) {
2566 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
2567
2568 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
2569
2570 (*this->pointsRep_)[counter][2] = coorRec[2] + t*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
2571
2572 pointsRepGlobMapping[counter] = r
2573 + s*(nmbPoints_oneDir_allSubdomain);
2574 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2575 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
2576 }
2577 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
2578 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
2579 }
2580
2581 pointsRepGlobMapping[counter] += t*nmbPoints_oneDir_allSubdomain*nmbPoints_oneDir;
2582 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
2583 pointsRepGlobMapping[counter] -= t*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2584 }
2585
2586 pointsRepGlobMapping[counter] += offset_x*(M)
2587 + offset_y*( nmbPoints_oneDir_allSubdomain * M );
2588 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2589 pointsRepGlobMapping[counter] -= offset_y*M*(nmbPoints_oneDir-1);
2590 }
2591 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
2592 pointsRepGlobMapping[counter] -= offset_y*M*(nmbPoints_oneDir-1);
2593 }
2594
2595 pointsRepGlobMapping[counter] += offset_z * M * nmbPoints_oneDir_allSubdomain * nmbPoints_oneDir;
2596 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
2597 pointsRepGlobMapping[counter] -= offset_z*M*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2598 }
2599
2600 pointsRepGlobMapping[counter] += offset_Squares_x * M * nmbSubdomainsSquares_OneDir;
2601 if (offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2602 pointsRepGlobMapping[counter] -= M * nmbSubdomainsSquares_OneDir;
2603 }
2604 else if(offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=M){
2605 pointsRepGlobMapping[counter] -= M * nmbSubdomainsSquares_OneDir;
2606 }
2607
2608 pointsRepGlobMapping[counter] += offset_Squares_z * nmbPoints_oneDir_allSubdomain * ((M) * nmbSubdomainsSquares_OneDir+1) * M * nmbSubdomainsSquares_OneDir;
2609 if (offset_Squares_z > 0 ) {
2610 pointsRepGlobMapping[counter] -= (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2611 }
2612 counter++;
2613 }
2614 }
2615 }
2616
2617
2618 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
2619
2620 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
2621
2622 if (verbose) {
2623 std::cout << "-- Building P1 Unique Points ... " << std::flush;
2624 }
2625
2626 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(), std::vector<double>(3,0.0)));
2627 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
2628 LO index;
2629
2630 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
2631 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
2632 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
2633 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
2634 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
2635 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
2636 }
2637
2638 if (verbose) {
2639 std::cout << " done! --" << std::endl;
2640 }
2641
2642
2643 if (verbose) {
2644 std::cout << "-- Building P1 Elements ... " << std::flush;
2645 }
2646 counter = 0;
2647 for (int t=0; t < M; t++) {
2648 for (int s=0; s < M; s++) {
2649 for (int r=0; r < M; r++) {
2650 (*elementsVec)[counter][0] = r+1 + (M+1) * s + (M+1)*(M+1) * t ;
2651 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
2652 (*elementsVec)[counter][2] = r+1 + (M+1) * s + (M+1)*(M+1) * (t+1) ;
2653 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2654 counter++;
2655 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * (t+1) ;
2656 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
2657 (*elementsVec)[counter][2] = r+1 + (M+1) * s + (M+1)*(M+1) * (t+1) ;
2658 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2659 counter++;
2660 (*elementsVec)[counter][0] = r+1 + (M+1) * s + (M+1)*(M+1) * t ;
2661 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * t ;
2662 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
2663 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2664 counter++;
2665 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
2666 (*elementsVec)[counter][1] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
2667 (*elementsVec)[counter][2] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * t ;
2668 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2669 counter++;
2670 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
2671 (*elementsVec)[counter][1] = r + (M+1) * (s+1) + (M+1)*(M+1) * t ;
2672 (*elementsVec)[counter][2] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2673 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2674 counter++;
2675 (*elementsVec)[counter][0] = r + (M+1) * s + (M+1)*(M+1) * t ;
2676 (*elementsVec)[counter][1] = r + (M+1) * s + (M+1)*(M+1) * (t+1) ;
2677 (*elementsVec)[counter][2] = r + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2678 (*elementsVec)[counter][3] = r+1 + (M+1) * (s+1) + (M+1)*(M+1) * (t+1) ;
2679 counter++;
2680 }
2681 }
2682 }
2683 if (verbose) {
2684 std::cout << " done! --" << std::endl;
2685 }
2686 buildElementsClass(elementsVec);
2687 }
2688 else if(FEType == "P1-disc" || FEType == "P1-disc-global")
2689 buildP1_Disc_Q2_3DBFS( N, MM, numProcsCoarseSolve );
2690 else if(FEType == "P2"){
2691
2692 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
2693 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2694 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(10,-1)));
2695 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2696
2697 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2698
2699 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2700 int offset_Squares_y = 0;
2701 int offset_Squares_z = ((whichSquareSet+1) % 2);
2702
2703 int counter = 0;
2704 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2705 int offset_y = 0;
2706 int offset_z = 0;
2707 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N) {
2708 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2709 }
2710
2711 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))>=N*N ) {
2712 offset_z = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))/(N*(N));
2713 }
2714
2715 bool p1point;
2716 int p1_s = 0;
2717 int p1_r = 0;
2718 int p1_t = 0;
2719 for (int t=0; t < 2*(M+1)-1; t++) {
2720 for (int s=0; s < 2*(M+1)-1; s++) {
2721 for (int r=0; r < 2*(M+1)-1; r++) {
2722 p1point = false;
2723 if (s%2==0 && r%2==0 && t%2==0) {
2724 p1point = true;
2725 p1_s = s/2;
2726 p1_r = r/2;
2727 p1_t = t/2;
2728 }
2729 (*this->pointsRep_)[counter][0] = coorRec[0] + r*h/2. + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
2730
2731 (*this->pointsRep_)[counter][1] = coorRec[1] + s*h/2. + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
2732
2733 (*this->pointsRep_)[counter][2] = coorRec[2] + t*h/2. + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
2734
2735 pointsRepGlobMapping[counter] = r
2736 + s*(nmbPoints_oneDir_allSubdomain);
2737 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2738 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
2739 }
2740 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
2741 pointsRepGlobMapping[counter] -= s*(nmbPoints_oneDir-1);
2742 }
2743
2744 pointsRepGlobMapping[counter] += t*nmbPoints_oneDir_allSubdomain*nmbPoints_oneDir;
2745 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
2746 pointsRepGlobMapping[counter] -= t*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2747 }
2748
2749 pointsRepGlobMapping[counter] += offset_x*(2*M)
2750 + offset_y*( nmbPoints_oneDir_allSubdomain * 2*M );
2751 if (offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2752 pointsRepGlobMapping[counter] -= offset_y*2*M*(nmbPoints_oneDir-1);
2753 }
2754 else if(offset_Squares_x > 0 && offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
2755 pointsRepGlobMapping[counter] -= offset_y*2*M*(nmbPoints_oneDir-1);
2756 }
2757
2758 pointsRepGlobMapping[counter] += offset_z * 2*M * nmbPoints_oneDir_allSubdomain * nmbPoints_oneDir;
2759 if (offset_Squares_x > 0 && offset_Squares_z == 0 ) {
2760 pointsRepGlobMapping[counter] -= offset_z*2*M*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2761 }
2762
2763 pointsRepGlobMapping[counter] += offset_Squares_x * 2*M * nmbSubdomainsSquares_OneDir;
2764 if (offset_Squares_z == 0 && offset_z+1!=nmbSubdomainsSquares_OneDir) {
2765 pointsRepGlobMapping[counter] -= 2*M * nmbSubdomainsSquares_OneDir;
2766 }
2767 else if(offset_Squares_z == 0 && offset_z+1==nmbSubdomainsSquares_OneDir && t!=2*M){
2768 pointsRepGlobMapping[counter] -= 2*M * nmbSubdomainsSquares_OneDir;
2769 }
2770
2771 pointsRepGlobMapping[counter] += offset_Squares_z * nmbPoints_oneDir_allSubdomain * ((2*M) * nmbSubdomainsSquares_OneDir+1) * 2*M * nmbSubdomainsSquares_OneDir;
2772 if (offset_Squares_z > 0 ) {
2773 pointsRepGlobMapping[counter] -= (nmbPoints_oneDir-1)*(nmbPoints_oneDir-1)*(nmbPoints_oneDir);
2774 }
2775 counter++;
2776 }
2777 }
2778 }
2779
2780 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
2781
2782 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
2783
2784 if (verbose) {
2785 std::cout << "-- Building P2 Unique Points ... " << std::flush;
2786 }
2787
2788 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
2789 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),0));
2790
2791 LO index;
2792 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
2793
2794 index = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement(i) );
2795
2796 (*this->pointsUni_)[i][0] = (*this->pointsRep_)[index][0];
2797 (*this->pointsUni_)[i][1] = (*this->pointsRep_)[index][1];
2798 (*this->pointsUni_)[i][2] = (*this->pointsRep_)[index][2];
2799 (*this->bcFlagUni_)[i] = (*this->bcFlagRep_)[index];
2800 }
2801
2802 if (verbose) {
2803 std::cout << " done! --" << std::endl;
2804 }
2805
2806 // Face 1 Face2 Face 3 Face 4
2807 // 2 2 * * 9 * * 3 3 * * 9 * * 2 3
2808 // * * * * * * * *
2809 // * * * * * * * *
2810 // 5 6 6 7 8 5 8 7
2811 // * * * * * * * *
2812 // * * * * * * * *
2813 // 1 * * 4 * * 0 0 1 1 * * 4 * * 0
2814
2815
2816 int P2M = 2*(M+1)-1;
2817
2818 counter = 0;
2819 for (int t=0; t < M; t++) {
2820 for (int s=0; s < M; s++) {
2821 for (int r=0; r < M; r++) {
2822
2823 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2824 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2825 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
2826 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2827
2828 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2829 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M;
2830 (*elementsVec)[counter][7] = 2*(r+1) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M;
2831 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M;
2832 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M;
2833 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1);
2834
2835 counter++;
2836
2837 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
2838 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2839 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
2840 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2841
2842 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
2843 (*elementsVec)[counter][6] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
2844 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t+1) ;
2845 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
2846 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2847 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t+1);
2848
2849 counter++;
2850
2851 (*elementsVec)[counter][0] = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2852 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2853 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
2854 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2855
2856 (*elementsVec)[counter][4] = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2857 (*elementsVec)[counter][6] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) ;
2858 (*elementsVec)[counter][7] = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2859 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) ;
2860 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2861 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
2862
2863 counter++;
2864
2865 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2866 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
2867 (*elementsVec)[counter][2] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
2868 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2869
2870 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
2871 (*elementsVec)[counter][6] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
2872 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2873 (*elementsVec)[counter][5] = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t);
2874 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
2875 (*elementsVec)[counter][9] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
2876
2877 counter++;
2878
2879 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2880 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) ;
2881 (*elementsVec)[counter][2] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2882 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2883
2884 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) ;
2885 (*elementsVec)[counter][6] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2886 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2887 (*elementsVec)[counter][5] = 2*(r) + 2*P2M*(s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
2888 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t) +P2M*P2M ;
2889 (*elementsVec)[counter][9] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t+1) ;
2890
2891 counter++;
2892
2893 (*elementsVec)[counter][0] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) ;
2894 (*elementsVec)[counter][1] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) ;
2895 (*elementsVec)[counter][2] = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2896 (*elementsVec)[counter][3] = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) ;
2897
2898 (*elementsVec)[counter][4] = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) +P2M*P2M ;
2899 (*elementsVec)[counter][6] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2900 (*elementsVec)[counter][7] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M ;
2901 (*elementsVec)[counter][5] = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) ;
2902 (*elementsVec)[counter][8] = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) ;
2903 (*elementsVec)[counter][9] = 2*(r) +1 + 2*P2M*(s+1) + 2*P2M*P2M * (t+1) ;
2904
2905 counter++;
2906
2907 }
2908 }
2909 }
2910 buildElementsClass(elementsVec);
2911 }
2912 else if(FEType == "Q1")
2913 build3DQ1BFS( N, MM, numProcsCoarseSolve);
2914 else if(FEType == "Q2")
2915 build3DQ2BFS( N, MM, numProcsCoarseSolve);
2916
2917};
2918
2919template <class SC, class LO, class GO, class NO>
2920void MeshStructured<SC,LO,GO,NO>::buildP1_Disc_Q2_3DBFS(int N,
2921 int M,
2922 int numProcsCoarseSolve){
2923
2924
2925
2926
2927 using Teuchos::RCP;
2928 using Teuchos::rcp;
2929 using Teuchos::ScalarTraits;
2930 typedef ScalarTraits<SC> ST;
2931
2932 bool verbose (this->comm_->getRank() == 0);
2933
2934 setRankRange( numProcsCoarseSolve );
2935
2936 if (verbose)
2937 std::cout << std::endl;
2938
2939 SC eps = ScalarTraits<SC>::eps();
2940
2941 int rank = this->comm_->getRank();
2942 int size = this->comm_->getSize() - numProcsCoarseSolve;
2943
2944
2945 int bfs_multiplier = (int) 2*(length)-1;
2946
2947 int nmbSubdomainsSquares = size / bfs_multiplier;
2948 int nmbSubdomainsSquares_OneDir = (std::pow(nmbSubdomainsSquares,1./3.) + 100*eps); // same as N
2949
2950 SC h = ST::one()/(M*N);
2951 SC H = ST::one()/N;
2952
2953 LO nmbElements = M*M*M;
2954 LO nmbPoints = 4*M*M*M; // 4 points for each element
2955
2956
2957 if (rank>=size) {
2958 M = -1; // keine Schleife wird ausgefuehrt
2959 nmbElements = 0;
2960 nmbPoints = 0;
2961 }
2962
2963 this->numElementsGlob_ = nmbElements * size;
2964
2965 int whichSquareSet = (int)rank / nmbSubdomainsSquares;
2966
2967 int offset_Squares_x = (int) (whichSquareSet+1) / 2;
2968 int offset_Squares_y = 0;
2969 int offset_Squares_z = ((whichSquareSet+1) % 2);
2970
2971 int counter = 0;
2972 int offset_x = ((rank - nmbSubdomainsSquares*whichSquareSet) % N);
2973 int offset_y = 0;
2974 int offset_z = 0;
2975 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))>=N)
2976 offset_y = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N))/(N);
2977
2978 if (((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))>=N*N )
2979 offset_z = (int) ((rank - nmbSubdomainsSquares*whichSquareSet) % (N*N*N))/(N*(N));
2980
2981
2982 this->pointsRep_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
2983 this->pointsUni_.reset(new std::vector<std::vector<double> >(nmbPoints,std::vector<double>(3,0.0)));
2984 this->bcFlagRep_.reset(new std::vector<int> (nmbPoints,0));
2985 this->bcFlagUni_.reset(new std::vector<int> (nmbPoints,0));
2986 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
2987
2988 if (verbose)
2989 std::cout << "-- Building P1-disc Points and Elements according to Q2 ... " << std::flush;
2990
2991 vec2D_int_ptr_Type elementsVec = Teuchos::rcp(new std::vector<std::vector<int> >(nmbElements,std::vector<int>(4,-1)));
2992
2993 counter = 0;
2994 LO pCounter = 0;
2995 GO globalCounterPoints = rank * nmbPoints;
2996 for (int t=0; t < M; t++) {
2997 for (int s=0; s < M; s++) {
2998 for (int r=0; r < M; r++) {
2999
3000 //point 1
3001 (*this->pointsRep_)[pCounter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
3002 (*this->pointsRep_)[pCounter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
3003 (*this->pointsRep_)[pCounter][2] = coorRec[2] + t*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
3004
3005 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
3006 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
3007 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
3008
3009 (*this->bcFlagRep_)[counter] = 1;
3010 }
3011
3012 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
3013 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
3014 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
3015
3016 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
3017
3018 (*elementsVec)[counter][0] = pCounter;
3019 pointsRepGlobMapping[pCounter] = globalCounterPoints;
3020 globalCounterPoints++;
3021 pCounter++;
3022
3023 //point 2
3024 (*this->pointsRep_)[pCounter][0] = coorRec[0] + (r+1)*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
3025 (*this->pointsRep_)[pCounter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
3026 (*this->pointsRep_)[pCounter][2] = coorRec[2] + t*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
3027 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
3028 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
3029 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
3030
3031 (*this->bcFlagRep_)[counter] = 1;
3032 }
3033
3034 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
3035 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
3036 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
3037
3038 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
3039
3040 (*elementsVec)[counter][1] = pCounter;
3041 pointsRepGlobMapping[pCounter] = globalCounterPoints;
3042 globalCounterPoints++;
3043 pCounter++;
3044 //point 3
3045 (*this->pointsRep_)[pCounter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
3046 (*this->pointsRep_)[pCounter][1] = coorRec[1] + (s+1)*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
3047 (*this->pointsRep_)[pCounter][2] = coorRec[2] + t*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
3048 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
3049 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
3050 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
3051
3052 (*this->bcFlagRep_)[counter] = 1;
3053 }
3054
3055 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
3056 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
3057 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
3058
3059 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
3060
3061 (*elementsVec)[counter][2] = pCounter;
3062 pointsRepGlobMapping[pCounter] = globalCounterPoints;
3063 globalCounterPoints++;
3064 pCounter++;
3065
3066 //point 4
3067 (*this->pointsRep_)[pCounter][0] = coorRec[0] + r*h + offset_x * H + offset_Squares_x * H * nmbSubdomainsSquares_OneDir;
3068 (*this->pointsRep_)[pCounter][1] = coorRec[1] + s*h + offset_y * H + offset_Squares_y * H * nmbSubdomainsSquares_OneDir;
3069 (*this->pointsRep_)[pCounter][2] = coorRec[2] + (t+1)*h + offset_z * H + offset_Squares_z * H * nmbSubdomainsSquares_OneDir;
3070 if ((*this->pointsRep_)[pCounter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[pCounter][0] < (coorRec[0]+eps) ||
3071 (*this->pointsRep_)[pCounter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[pCounter][1] < (coorRec[1]+eps) ||
3072 (*this->pointsRep_)[pCounter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[pCounter][2] < (coorRec[2]+eps) ) {
3073
3074 (*this->bcFlagRep_)[counter] = 1;
3075 }
3076
3077 (*this->pointsUni_)[pCounter][0] = (*this->pointsRep_)[pCounter][0];
3078 (*this->pointsUni_)[pCounter][1] = (*this->pointsRep_)[pCounter][1];
3079 (*this->pointsUni_)[pCounter][2] = (*this->pointsRep_)[pCounter][2];
3080
3081 (*this->bcFlagUni_)[pCounter] = (*this->bcFlagRep_)[pCounter];
3082
3083 (*elementsVec)[counter][3] = pCounter;
3084 pointsRepGlobMapping[pCounter] = globalCounterPoints;
3085 globalCounterPoints++;
3086 pCounter++;
3087
3088 counter++;
3089 }
3090 }
3091 }
3092 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
3093
3094 this->mapUnique_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
3095
3096 buildElementsClass(elementsVec);
3097
3098 if (verbose)
3099 std::cout << "done!" << std::endl;
3100
3101}
3102template <class SC, class LO, class GO, class NO>
3103void MeshStructured<SC,LO,GO,NO>::setStructuredMeshFlags(int flagsOption,std::string FEType){
3104
3105 double tol=1.e-12;
3106
3107 switch (this->dim_) {
3108 case 2:
3109 switch (flagsOption) {
3110 case 0:
3111 break;
3112 case 1: //Rectangle left inflow
3113 for (int i=0; i<this->pointsUni_->size(); i++) {
3114 if (this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + height - tol)) {
3115 this->bcFlagUni_->at(i) = 3; //outflow
3116 }
3117 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol)) {
3118 this->bcFlagUni_->at(i) = 2; //inflow
3119 }
3120 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3121 this->bcFlagUni_->at(i) = 1;
3122 }
3123 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3124 this->bcFlagUni_->at(i) = 1;
3125 }
3126 }
3127 for (int i=0; i<this->pointsRep_->size(); i++) {
3128 if (this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + height - tol)) {
3129 this->bcFlagRep_->at(i) = 3;
3130 }
3131 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol)) {
3132 this->bcFlagRep_->at(i) = 2;
3133 }
3134 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3135 this->bcFlagRep_->at(i) = 1;
3136 }
3137 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3138 this->bcFlagRep_->at(i) = 1;
3139 }
3140 }
3141 break;
3142 case 2: //BFS
3143 for (int i=0; i<this->pointsUni_->size(); i++) {
3144 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol) && this->pointsUni_->at(i).at(1) > (coorRec[1]+1. -tol) && this->pointsUni_->at(i).at(1) < (coorRec[1]+ height +tol)) {
3145 this->bcFlagUni_->at(i) = 2;
3146 }
3147 if (this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) || this->pointsUni_->at(i).at(1) > (coorRec[1]+height - tol) || this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol)) {
3148 this->bcFlagUni_->at(i) = 1;
3149 }
3150 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(0) < (coorRec[0]+1. + tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + 1. - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + 1. + tol)) {
3151 this->bcFlagUni_->at(i) = 1;
3152 }
3153 if (this->pointsUni_->at(i).at(1) > (coorRec[1] - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + 1. + tol) && this->pointsUni_->at(i).at(0) > (coorRec[0] + 1. - tol) && this->pointsUni_->at(i).at(0) < (coorRec[0] + 1. + tol)) {
3154 this->bcFlagUni_->at(i) = 1;
3155 }
3156 if (this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + height - tol)) {
3157 this->bcFlagUni_->at(i) = 3;
3158 }
3159 }
3160 for (int i=0; i<this->pointsRep_->size(); i++) {
3161 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol) && this->pointsRep_->at(i).at(1) > (coorRec[1]+1. -tol) && this->pointsRep_->at(i).at(1) < (coorRec[1]+ height +tol)) {
3162 this->bcFlagRep_->at(i) = 2;
3163 }
3164 if (this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) || this->pointsRep_->at(i).at(1) > (coorRec[1]+height - tol) || this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol)) {
3165 this->bcFlagRep_->at(i) = 1;
3166 }
3167 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(0) < (coorRec[0]+1. + tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + 1. - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + 1. + tol)) {
3168 this->bcFlagRep_->at(i) = 1;
3169 }
3170 if (this->pointsRep_->at(i).at(1) > (coorRec[1] - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + 1. + tol) && this->pointsRep_->at(i).at(0) > (coorRec[0] + 1. - tol) && this->pointsRep_->at(i).at(0) < (coorRec[0] + 1. + tol)) {
3171 this->bcFlagRep_->at(i) = 1;
3172 }
3173 if (this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + height - tol)) {
3174 this->bcFlagRep_->at(i) = 3;
3175 }
3176 }
3177 break;
3178 case 3: //tpm
3179 for (int i=0; i<this->pointsUni_->size(); i++) {
3180 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol)) {
3181 this->bcFlagUni_->at(i) = 2; //left
3182 }
3183 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3184 this->bcFlagUni_->at(i) = 1; //bottom
3185 }
3186 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3187 this->bcFlagUni_->at(i) = 4; //top
3188 }
3189 if (this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + height - tol)) {
3190 this->bcFlagUni_->at(i) = 3; //right
3191 }
3192 }
3193 for (int i=0; i<this->pointsRep_->size(); i++) {
3194 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol)) {
3195 this->bcFlagRep_->at(i) = 2;
3196 }
3197 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3198 this->bcFlagRep_->at(i) = 1;
3199 }
3200 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3201 this->bcFlagRep_->at(i) = 4;
3202 }
3203 if (this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + height - tol)) {
3204 this->bcFlagRep_->at(i) = 3;
3205 }
3206 }
3207 break;
3208 case 4: //mini tpm, we set all values manually
3209 if (FEType == "P2") {
3210 this->bcFlagUni_->at(0) = 1;
3211 this->bcFlagUni_->at(1) = 2;
3212 this->bcFlagUni_->at(2) = 2;
3213 this->bcFlagUni_->at(3) = 2;
3214 this->bcFlagUni_->at(4) = 4; // Dirichlet and lineload
3215 this->bcFlagUni_->at(5) = 3;
3216 this->bcFlagUni_->at(6) = 0;
3217 this->bcFlagUni_->at(7) = 0;
3218 this->bcFlagUni_->at(8) = 0;
3219 this->bcFlagUni_->at(9) = 5; // lineload
3220 this->bcFlagUni_->at(10) = 1;
3221 this->bcFlagUni_->at(11) = 2;
3222 this->bcFlagUni_->at(12) = 2;
3223 this->bcFlagUni_->at(13) = 2;
3224 this->bcFlagUni_->at(14) = 4;
3225
3226 this->bcFlagRep_->at(0) = 1;
3227 this->bcFlagRep_->at(1) = 2;
3228 this->bcFlagRep_->at(2) = 2;
3229 this->bcFlagRep_->at(3) = 2;
3230 this->bcFlagRep_->at(4) = 4;
3231 this->bcFlagRep_->at(5) = 3;
3232 this->bcFlagRep_->at(6) = 0;
3233 this->bcFlagRep_->at(7) = 0;
3234 this->bcFlagRep_->at(8) = 0;
3235 this->bcFlagRep_->at(9) = 5;
3236 this->bcFlagRep_->at(10) = 1;
3237 this->bcFlagRep_->at(11) = 2;
3238 this->bcFlagRep_->at(12) = 2;
3239 this->bcFlagRep_->at(13) = 2;
3240 this->bcFlagRep_->at(14) = 4;
3241 } else if(FEType=="P1") {
3242 this->bcFlagUni_->at(0) = 1;
3243 this->bcFlagUni_->at(1) = 1;
3244 this->bcFlagUni_->at(2) = 1;
3245 this->bcFlagUni_->at(3) = 2;
3246 this->bcFlagUni_->at(4) = 2;
3247 this->bcFlagUni_->at(5) = 1;
3248
3249 this->bcFlagRep_->at(0) = 1;
3250 this->bcFlagRep_->at(1) = 1;
3251 this->bcFlagRep_->at(2) = 1;
3252 this->bcFlagRep_->at(3) = 2;
3253 this->bcFlagRep_->at(4) = 2;
3254 this->bcFlagRep_->at(5) = 1;
3255 }
3256 break;
3257 case 5: // LDC Square
3258 for (int i=0; i<this->pointsUni_->size(); i++) {
3259 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3260 this->bcFlagUni_->at(i) = 1; // bottom
3261 }
3262 if (this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + height - tol)) {
3263 this->bcFlagUni_->at(i) = 1; // right
3264 }
3265 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol)) {
3266 this->bcFlagUni_->at(i) = 1; //Left
3267 }
3268 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3269 this->bcFlagUni_->at(i) = 2; // top
3270 }
3271 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] +tol)) {
3272 this->bcFlagUni_->at(i) = 3; // (0,0) point of ldc
3273 }
3274 }
3275 for (int i=0; i<this->pointsRep_->size(); i++) {
3276 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3277 this->bcFlagRep_->at(i) = 1; // bottom
3278 }
3279 if (this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + height - tol)) {
3280 this->bcFlagRep_->at(i) = 1;
3281 }
3282 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol)) {
3283 this->bcFlagRep_->at(i) = 1;
3284 }
3285 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + height - tol) ) {
3286 this->bcFlagRep_->at(i) = 2;
3287 }
3288 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] +tol)) {
3289 this->bcFlagRep_->at(i) = 3;
3290 }
3291 }
3292 break;
3293 default:
3294 break;
3295 }
3296 break;
3297 case 3:
3298 switch (flagsOption) {
3299 case 0:
3300 break;
3301 case 1: //Rectangle left inflow
3302 for (int i=0; i<this->pointsUni_->size(); i++) {
3303 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) ) {
3304 this->bcFlagUni_->at(i) = 2;
3305 }
3306 //out
3307 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + length - tol) &&
3308 this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) &&
3309 this->pointsUni_->at(i).at(1) < (coorRec[1] + width - tol)&&
3310 this->pointsUni_->at(i).at(2) > (coorRec[2] + tol) &&
3311 this->pointsUni_->at(i).at(2) < (coorRec[2] + height - tol)) {
3312 this->bcFlagUni_->at(i) = 3;
3313 }
3314 //bottom
3315 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3316 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) ) {
3317 this->bcFlagUni_->at(i) = 1;
3318 }
3319 //top
3320 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3321 this->pointsUni_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3322 this->bcFlagUni_->at(i) = 1;
3323 }
3324 //front
3325 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3326 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3327 this->bcFlagUni_->at(i) = 1;
3328 }
3329 //back
3330 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3331 this->pointsUni_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3332 this->bcFlagUni_->at(i) = 1;
3333 }
3334 }
3335 for (int i=0; i<this->pointsUni_->size(); i++) {
3336 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) ) {
3337 this->bcFlagRep_->at(i) = 2;
3338 }
3339 //out
3340 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + length - tol) &&
3341 this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) &&
3342 this->pointsRep_->at(i).at(1) < (coorRec[1] + width - tol)&&
3343 this->pointsRep_->at(i).at(2) > (coorRec[2] + tol) &&
3344 this->pointsRep_->at(i).at(2) < (coorRec[2] + height - tol)) {
3345 this->bcFlagRep_->at(i) = 3;
3346 }
3347 //bottom
3348 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3349 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) ) {
3350 this->bcFlagRep_->at(i) = 1;
3351 }
3352 //top
3353 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3354 this->pointsRep_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3355 this->bcFlagRep_->at(i) = 1;
3356 }
3357 //front
3358 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3359 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3360 this->bcFlagRep_->at(i) = 1;
3361 }
3362 //back
3363 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3364 this->pointsRep_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3365 this->bcFlagRep_->at(i) = 1;
3366 }
3367 }
3368 break;
3369 case 2: //BFS
3370 for (int i=0; i<this->pointsUni_->size(); i++) {
3371 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] - tol) && this->pointsUni_->at(i).at(1) < (coorRec[1]+ width +tol)
3372 && this->pointsUni_->at(i).at(2) > (coorRec[2]+1. -tol) && this->pointsUni_->at(i).at(1) < (coorRec[2]+ height +tol)) {
3373 this->bcFlagUni_->at(i) = 2;
3374 }
3375 //bottom top
3376 if (this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) || (this->pointsUni_->at(i).at(2) < (coorRec[2]+1. + tol) && this->pointsUni_->at(i).at(0) < (coorRec[0]+1. + tol)) || this->pointsUni_->at(i).at(2) > (coorRec[2]+height - tol)) {
3377 this->bcFlagUni_->at(i) = 1;
3378 }
3379 //front back
3380 if (this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) || this->pointsUni_->at(i).at(1) > (coorRec[1]+width - tol)) {
3381 this->bcFlagUni_->at(i) = 1;
3382 }
3383 //step left
3384 if (this->pointsUni_->at(i).at(0) < (coorRec[0]+1. + tol) && this->pointsUni_->at(i).at(0) > (coorRec[0]+1. - tol) && this->pointsUni_->at(i).at(2) < (coorRec[2]+1.+tol)) {
3385 this->bcFlagUni_->at(i) = 1;
3386 }
3387 if (this->pointsUni_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] + width - tol)
3388 && this->pointsUni_->at(i).at(2) > (coorRec[2] + tol) && this->pointsUni_->at(i).at(2) < (coorRec[2] + height - tol)) {
3389 this->bcFlagUni_->at(i) = 3;
3390 }
3391 }
3392 for (int i=0; i<this->pointsRep_->size(); i++) {
3393 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] - tol) && this->pointsRep_->at(i).at(1) < (coorRec[1]+ width +tol)
3394 && this->pointsRep_->at(i).at(2) > (coorRec[2]+1. -tol) && this->pointsRep_->at(i).at(1) < (coorRec[2]+ height +tol)) {
3395 this->bcFlagRep_->at(i) = 2;
3396 }
3397 //bottom top
3398 if (this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) || (this->pointsRep_->at(i).at(2) < (coorRec[2]+1. + tol) && this->pointsRep_->at(i).at(0) < (coorRec[0]+1. + tol)) || this->pointsRep_->at(i).at(2) > (coorRec[2]+height - tol)) {
3399 this->bcFlagRep_->at(i) = 1;
3400 }
3401 //front back
3402 if (this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) || this->pointsRep_->at(i).at(1) > (coorRec[1]+width - tol)) {
3403 this->bcFlagRep_->at(i) = 1;
3404 }
3405 //step left
3406 if (this->pointsRep_->at(i).at(0) < (coorRec[0]+1. + tol) && this->pointsRep_->at(i).at(0) > (coorRec[0]+1. - tol) && this->pointsRep_->at(i).at(2) < (coorRec[2]+1.+tol)) {
3407 this->bcFlagRep_->at(i) = 1;
3408 }
3409 if (this->pointsRep_->at(i).at(0) > (coorRec[0]+length - tol) && this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] + width - tol)
3410 && this->pointsRep_->at(i).at(2) > (coorRec[2] + tol) && this->pointsRep_->at(i).at(2) < (coorRec[2] + height - tol)) {
3411 this->bcFlagRep_->at(i) = 3;
3412 }
3413 }
3414 break;
3415 case 3: //SMC Cube Test
3416 for (int i=0; i<this->pointsUni_->size(); i++) {
3417 // z=1 Face
3418 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3419 this->pointsUni_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3420 this->bcFlagUni_->at(i) = 4;
3421 }
3422 // y=1 Face
3423 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3424 this->pointsUni_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3425 this->bcFlagUni_->at(i) = 5;
3426 }
3427 // x=1 Face
3428 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + length - tol) &&
3429 this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) &&
3430 this->pointsUni_->at(i).at(1) < (coorRec[1] + width - tol)&&
3431 this->pointsUni_->at(i).at(2) > (coorRec[2] + tol) &&
3432 this->pointsUni_->at(i).at(2) < (coorRec[2] + height - tol)) {
3433 this->bcFlagUni_->at(i) = 6;
3434 }
3435 // x=0 Face
3436 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) ) {
3437 this->bcFlagUni_->at(i) = 1;
3438 }
3439 // y=0 Face
3440 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3441 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3442 this->bcFlagUni_->at(i) = 2;
3443 }
3444 // z=0 Face
3445 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3446 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) ) {
3447 this->bcFlagUni_->at(i) = 3;
3448 }
3449
3450 // (0,0,z) Edge
3451 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) &&
3452 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) )
3453 this->bcFlagUni_->at(i) = 7;
3454
3455 // (0,y,0) Edge
3456 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) &&
3457 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) )
3458 this->bcFlagUni_->at(i) = 9;
3459
3460 // (x,0,0) Edge
3461 if (this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) &&
3462 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) )
3463 this->bcFlagUni_->at(i) = 8;
3464
3465 // (0,0,0) Point
3466 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) &&
3467 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) &&
3468 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) )
3469 this->bcFlagUni_->at(i) = 0;
3470
3471 }
3472 for (int i=0; i<this->pointsRep_->size(); i++) {
3473
3474 // z=1 Face
3475 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3476 this->pointsRep_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3477 this->bcFlagRep_->at(i) = 4;
3478 }
3479
3480 // y=1 face
3481 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3482 this->pointsRep_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3483 this->bcFlagRep_->at(i) = 5;
3484 }
3485 //x=1
3486 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + length - tol) &&
3487 this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) &&
3488 this->pointsRep_->at(i).at(1) < (coorRec[1] + width - tol)&&
3489 this->pointsRep_->at(i).at(2) > (coorRec[2] + tol) &&
3490 this->pointsRep_->at(i).at(2) < (coorRec[2] + height - tol)) {
3491 this->bcFlagRep_->at(i) = 6;
3492 }
3493 // x=0 Face
3494 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) ) {
3495 this->bcFlagRep_->at(i) = 1;
3496 }
3497 // z=0 Face
3498 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3499 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) ) {
3500 this->bcFlagRep_->at(i) = 3;
3501 }
3502 // y=0 Face
3503 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3504 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3505 this->bcFlagRep_->at(i) = 2;
3506 }
3507 // (0,0,z) Edge
3508 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) &&
3509 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) )
3510 this->bcFlagRep_->at(i) = 7;
3511
3512 // (0,y,0) Edge
3513 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) &&
3514 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) )
3515 this->bcFlagRep_->at(i) = 9;
3516
3517 // (x,0,0) Edge
3518 if (this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) &&
3519 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) )
3520 this->bcFlagRep_->at(i) = 8;
3521 // (0,0,0) Point
3522 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) &&
3523 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) &&
3524 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) )
3525 this->bcFlagRep_->at(i) = 0;
3526
3527 }
3528 break;
3529 case 4: // tube flow through z-direction
3530 for (int i=0; i<this->pointsUni_->size(); i++) {
3531 //bottom
3532 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3533 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) ) {
3534 this->bcFlagUni_->at(i) = 4;
3535 }
3536 //top
3537 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3538 this->pointsUni_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3539 this->bcFlagUni_->at(i) = 5;
3540 }
3541 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) ) {
3542 this->bcFlagUni_->at(i) = 6;
3543 }
3544 //front
3545 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3546 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3547 this->bcFlagUni_->at(i) = 6;
3548 }
3549 //back
3550 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + tol) &&
3551 this->pointsUni_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3552 this->bcFlagUni_->at(i) = 6;
3553 }
3554 //out
3555 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + length - tol) &&
3556 this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) &&
3557 this->pointsUni_->at(i).at(1) < (coorRec[1] + width - tol)&&
3558 this->pointsUni_->at(i).at(2) > (coorRec[2] - tol) &&
3559 this->pointsUni_->at(i).at(2) < (coorRec[2] + height + tol)) {
3560 this->bcFlagUni_->at(i) = 6;
3561 }
3562 }
3563 for (int i=0; i<this->pointsUni_->size(); i++) {
3564
3565 //bottom
3566 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3567 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) ) {
3568 this->bcFlagRep_->at(i) = 4;
3569 }
3570 if (this->pointsRep_->at(i).at(0) < (coorRec[0] + tol) ) {
3571 this->bcFlagRep_->at(i) = 6;
3572 }
3573 //front
3574 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + tol) &&
3575 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3576 this->bcFlagRep_->at(i) = 6;
3577 }
3578 //back
3579 if (this->pointsRep_->at(i).at(0) >= (coorRec[0] + tol) &&
3580 this->pointsRep_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3581 this->bcFlagRep_->at(i) = 6;
3582 }
3583 //out
3584 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + length - tol) &&
3585 this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) &&
3586 this->pointsRep_->at(i).at(1) < (coorRec[1] + width - tol)&&
3587 this->pointsRep_->at(i).at(2) > (coorRec[2] - tol) &&
3588 this->pointsRep_->at(i).at(2) < (coorRec[2] + height + tol)) {
3589 this->bcFlagRep_->at(i) = 6;
3590 }
3591 //top
3592 if (this->pointsRep_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3593 this->bcFlagRep_->at(i) = 5;
3594 }
3595 }
3596 break;
3597 case 5: //LDC
3598 for (int i=0; i<this->pointsUni_->size(); i++) {
3599 if (this->pointsUni_->at(i).at(0) < (coorRec[0] + tol) ) {
3600 this->bcFlagUni_->at(i) = 1;
3601 }
3602 //bottom
3603 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3604 this->pointsUni_->at(i).at(2) < (coorRec[2] + tol) ) {
3605 this->bcFlagUni_->at(i) = 1;
3606 }
3607 //front
3608 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3609 this->pointsUni_->at(i).at(1) < (coorRec[1] + tol) ) {
3610 this->bcFlagUni_->at(i) = 1;
3611 }
3612 //back
3613 if (this->pointsUni_->at(i).at(0) > (coorRec[0] - tol) &&
3614 this->pointsUni_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3615 this->bcFlagUni_->at(i) = 1;
3616 }
3617 //out
3618 if (this->pointsUni_->at(i).at(0) > (coorRec[0] + length - tol) &&
3619 this->pointsUni_->at(i).at(1) > (coorRec[1] + tol) &&
3620 this->pointsUni_->at(i).at(1) < (coorRec[1] + width - tol)&&
3621 this->pointsUni_->at(i).at(2) > (coorRec[2] + tol) &&
3622 this->pointsUni_->at(i).at(2) < (coorRec[2] + height - tol)) {
3623 this->bcFlagUni_->at(i) = 1;
3624 }
3625 //top
3626 if (this->pointsUni_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3627 this->bcFlagUni_->at(i) = 2;
3628 }
3629 if (this->pointsUni_->at(i).at(0) < (coorRec[0] +tol) && this->pointsUni_->at(i).at(1) < (coorRec[1] +tol) && this->pointsUni_->at(i).at(2) < (coorRec[2] +tol)) {
3630 this->bcFlagUni_->at(i) = 3; // (0,0) point of ldc
3631 }
3632 }
3633 for (int i=0; i<this->pointsUni_->size(); i++) {
3634 if (this->pointsRep_->at(i).at(0) < (coorRec[0] - tol) ) {
3635 this->bcFlagRep_->at(i) = 1;
3636 }
3637 //bottom
3638 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3639 this->pointsRep_->at(i).at(2) < (coorRec[2] + tol) ) {
3640 this->bcFlagRep_->at(i) = 1;
3641 }
3642 //top
3643 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3644 this->pointsRep_->at(i).at(2) > (coorRec[2] + height - tol) ) {
3645 this->bcFlagRep_->at(i) = 2;
3646 }
3647 //front
3648 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3649 this->pointsRep_->at(i).at(1) < (coorRec[1] + tol) ) {
3650 this->bcFlagRep_->at(i) = 1;
3651 }
3652 //back
3653 if (this->pointsRep_->at(i).at(0) > (coorRec[0] - tol) &&
3654 this->pointsRep_->at(i).at(1) > (coorRec[1] + width - tol) ) {
3655 this->bcFlagRep_->at(i) = 1;
3656 }
3657 //out
3658 if (this->pointsRep_->at(i).at(0) > (coorRec[0] + length - tol) &&
3659 this->pointsRep_->at(i).at(1) > (coorRec[1] + tol) &&
3660 this->pointsRep_->at(i).at(1) < (coorRec[1] + width - tol)&&
3661 this->pointsRep_->at(i).at(2) > (coorRec[2] + tol) &&
3662 this->pointsRep_->at(i).at(2) < (coorRec[2] + height - tol)) {
3663 this->bcFlagRep_->at(i) = 1;
3664 }
3665 if (this->pointsRep_->at(i).at(0) < (coorRec[0] +tol) && this->pointsRep_->at(i).at(1) < (coorRec[1] +tol) && this->pointsRep_->at(i).at(2) < (coorRec[2] +tol)) {
3666 this->bcFlagRep_->at(i) = 3; // (0,0) point of ldc
3667 }
3668 }
3669 break;
3670 default:
3671 break;
3672 }
3673
3674 break;
3675 default:
3676 break;
3677 }
3678}
3679
3680template <class SC, class LO, class GO, class NO>
3681void MeshStructured<SC,LO,GO,NO>::buildMesh3D5Elements(std::string FEType,
3682 int N,
3683 int M,
3684 int numProcsCoarseSolve){
3685
3686 using Teuchos::RCP;
3687 using Teuchos::rcp;
3688 using Teuchos::ScalarTraits;
3689
3690 TEUCHOS_TEST_FOR_EXCEPTION(!(M>=1),std::logic_error,"H/h is to small.");
3691 TEUCHOS_TEST_FOR_EXCEPTION(this->comm_.is_null(),std::runtime_error,"comm_ is null.");
3692
3693 bool verbose (this->comm_->getRank() == 0);
3694
3695 if (verbose) {
3696 std::cout << "-- Building structured 3D Mesh --" << std::endl;
3697 }
3698
3699 setRankRange( numProcsCoarseSolve );
3700
3701 SC eps = ScalarTraits<SC>::eps();
3702
3703 int rank = this->comm_->getRank();
3704 int size = this->comm_->getSize() - numProcsCoarseSolve;
3705
3706 SC h = length/(M*N);//add variable length/width/heigth
3707 SC H = length/N;
3708
3709 if (verbose) {
3710 std::cout << "-- H:"<<H << " h:" <<h << " --" << std::endl;
3711 std::cout << "-- N:"<<N << " M:" <<M << " --" << std::endl;
3712 }
3713
3714 LO nmbPoints_oneDir;
3715
3716 LO nmbElements;
3717 LO nmbPoints;
3718 if (FEType == "P2"){
3719 nmbPoints_oneDir = N * (2*(M+1)-1) - (N-1);
3720 nmbPoints = (2*(M+1)-1)*(2*(M+1)-1)*(2*(M+1)-1) - M*M*M;
3721 }
3722 else
3723 TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error, "Wrong FE-Type, 5 element subcube devision only available with P2 discretization.");
3724
3725
3726
3727 if (verbose) {
3728 std::cout << "-- Number of Points in one direction: " << nmbPoints_oneDir << " || Number of Points " << nmbPoints << " --" << std::endl;
3729 }
3730 this->FEType_ = FEType;
3731
3732 GO nmbPGlob_oneDir = N * (M+1) - (N-1);
3733
3734 this->numElementsGlob_ = 5*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1)*(nmbPGlob_oneDir-1);
3735 int MM=M;
3736 if (rank>=size) {
3737 M = -1; // keine Schleife wird ausgefuehrt
3738 nmbElements = 0;
3739 nmbPoints = 0;
3740 }
3741 else{
3742 nmbElements = 5*M*M*M; //6*M*M*M;
3743 }
3744 vec2D_int_ptr_Type elementsVec;
3745 vec_int_ptr_Type elementFlag;
3746
3747
3748 this->pointsRep_.reset(new vec2D_dbl_Type(nmbPoints, vec_dbl_Type(3, 0.0)));
3749 this->bcFlagRep_.reset(new vec_int_Type (nmbPoints, 10));
3750 elementsVec = Teuchos::rcp(new vec2D_int_Type(nmbElements, vec_int_Type(10, -1)));
3751 elementFlag = Teuchos::rcp(new vec_int_Type( elementsVec->size(),0 ) );
3752 Teuchos::Array<GO> pointsRepGlobMapping(nmbPoints);
3753
3754 int counter = 0;
3755 int offset_x = (rank % N); // -M*M*M;
3756 int offset_y = 0;
3757 int offset_z = 0;
3758
3759 if ((rank % (N*N))>=N) {
3760 offset_y = (int) (rank % (N*N))/(N);
3761 }
3762
3763 if ((rank % (N*N*N))>=N*N ) {
3764 offset_z = (int) (rank % (N*N*N))/(N*(N));
3765 }
3766
3767 if (verbose) {
3768 std::cout << "-- Building P2 Points Repeated ... " << std::endl;
3769 }
3770 std::cout << " Offsets on Rank " << rank << " || x=" << offset_x << " y=" << offset_y << " z=" << offset_z << std::endl;
3771 this->comm_->barrier();
3772
3773 bool p1point;
3774 int p1_s = 0;
3775 int p1_r = 0;
3776 int p1_t = 0;
3777 int nodeSkip = 0;
3778 int offset=0;
3779 for (int t=0; t < 2*(M+1)-1; t++) {
3780 for (int s=0; s < 2*(M+1)-1; s++) {
3781 for (int r=0; r < 2*(M+1)-1; r++) {
3782
3783 p1point = false;
3784 if (s%2==0 && r%2==0 && t%2==0) {
3785 p1point = true;
3786 p1_s = s/2;
3787 p1_r = r/2;
3788 p1_t = t/2;
3789 }
3790
3791 (*this->pointsRep_)[counter][0] = r*h/2.0 + offset_x * H;
3792 if ((*this->pointsRep_)[counter][0]<eps && (*this->pointsRep_)[counter][0]>-eps) (*this->pointsRep_)[counter][0]=0.0;
3793 (*this->pointsRep_)[counter][1] = s*h/2.0 + offset_y * H;
3794 if ((*this->pointsRep_)[counter][1]<eps && (*this->pointsRep_)[counter][1]>-eps) (*this->pointsRep_)[counter][1]=0.0;
3795 (*this->pointsRep_)[counter][2] = t*h/2.0 + offset_z * H;
3796 if ((*this->pointsRep_)[counter][2]<eps && (*this->pointsRep_)[counter][2]>-eps) (*this->pointsRep_)[counter][2]=0.0;
3797
3798 if(N>1){
3799 offset=0;
3800 if(s%2 == 1 && t%2 == 1 && r%2==0)
3801 offset = M*offset_x +N*M*M*offset_y + N*M*(s-1)/2; // pre-skip
3802
3803 if(s%2 == 0 && t%2==1 ){ // letzte wand der subwürfel nach skip
3804 offset += M*M*N*offset_y + s/2*M*N;//M*offset_x
3805 nodeSkip=0;
3806 }
3807 if(t%2==0 && t>0 )
3808 offset += M*M*N*N * t/2; // z- ebenen zwischen subdomains
3809 if(t%2==1 && t>0 )
3810 offset += M*M*N*N * (t-1)/2;
3811 }
3812
3813 pointsRepGlobMapping[counter] = r + s*nmbPoints_oneDir + t*nmbPoints_oneDir*nmbPoints_oneDir \
3814 + offset_x*(2*(M+1)-2) + offset_y*(nmbPoints_oneDir)*(2*(M+1)-2) + offset_z*(nmbPoints_oneDir)*(nmbPoints_oneDir)*(2*(M+1)-2) \
3815 -offset_z*(N*N*M*M*M)-nodeSkip-offset;
3816
3817 if ((*this->pointsRep_)[counter][0] > (coorRec[0]+length-eps) || (*this->pointsRep_)[counter][0] < (coorRec[0]+eps) ||
3818 (*this->pointsRep_)[counter][1] > (coorRec[1]+width-eps) || (*this->pointsRep_)[counter][1] < (coorRec[1]+eps) ||
3819 (*this->pointsRep_)[counter][2] > (coorRec[2]+height-eps) || (*this->pointsRep_)[counter][2] < (coorRec[2]+eps) ) {
3820 (*this->bcFlagRep_)[counter] = 1;
3821
3822 }
3823 if (s%2==1 && r%2==1 && t%2==1) {
3824 nodeSkip ++;
3825
3826 }
3827 else
3828 counter++;
3829
3830
3831 }
3832 }
3833 }
3834
3835 counter =0;
3836 int P2M = 2*(M+1)-1;
3837
3838 this->mapRepeated_.reset(new Map<LO,GO,NO>( (GO) -1, pointsRepGlobMapping(), 0, this->comm_) );
3839
3840
3841 this->mapUnique_ = this->mapRepeated_->buildUniqueMap( numProcsCoarseSolve );
3842
3843 this->pointsUni_.reset(new std::vector<std::vector<double> >(this->mapUnique_->getNodeNumElements(),std::vector<double>(3,0.0)));
3844 this->bcFlagUni_.reset(new std::vector<int> (this->mapUnique_->getNodeNumElements(),10));
3845
3846 if (verbose) {
3847 std::cout << "-- Building P2 Points Unique ... " << std::endl;
3848 }
3849 if (verbose) {
3850 std::cout << "-- Number of repeated points per proc: " << this->mapRepeated_->getNodeNumElements() << " ... " << std::endl;
3851 }
3852
3853 this->comm_->barrier();
3854
3855 // Points and Flags Unique
3856 this->pointsUni_.reset(new std::vector<std::vector<double> >( this->mapUnique_->getNodeNumElements(), std::vector<double>(this->dim_,-1. ) ) );
3857 this->bcFlagUni_.reset( new std::vector<int> ( this->mapUnique_->getNodeNumElements(), 0 ) );
3858 for (int i=0; i<this->mapUnique_->getNodeNumElements(); i++) {
3859 GO gid = this->mapUnique_->getGlobalElement( i );
3860 LO id = this->mapRepeated_->getLocalElement( this->mapUnique_->getGlobalElement( i ) );
3861 this->pointsUni_->at(i) = this->pointsRep_->at(id);
3862 this->bcFlagUni_->at(i) = this->bcFlagRep_->at(id);
3863
3864 }
3865 this->comm_->barrier();
3866
3867
3868
3869
3870 // Face 1 Face2 Face 3 Face 4
3871 // 2 2 * * 9 * * 3 3 * * 9 * * 2 3
3872 // * * * * * * * *
3873 // * * * * * * * *
3874 // 5 6 6 7 8 5 8 7
3875 // * * * * * * * *
3876 // * * * * * * * *
3877 // 1 * * 4 * * 0 0 1 1 * * 4 * * 0
3878
3879
3880 //int P2M = 2*(M+1)-1;
3881
3882 if (verbose) {
3883 std::cout << "-- ElementsList ... " << std::endl;
3884 std::cout << "-- P2M =" << P2M << std::endl;
3885 }
3886
3887
3888 counter = 0;
3889 nodeSkip=0;
3890 for (int t=0; t < M; t++) { // z
3891 for (int s=0; s < M; s++) { // y
3892 for (int r=0; r < M; r++) { // x
3893
3894
3895 // x-y = z = 0
3896 int n_0 = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t) -t*M*M;
3897 int n_1 = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t) -t*M*M;
3898 int n_2 = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t) -t*M*M ;
3899
3900 int n_3 = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) -t*M*M;
3901 int n_4 = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) -t*M*M;
3902 int n_5 = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) -t*M*M ;
3903
3904 int n_6 = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) -t*M*M;
3905 int n_7 = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) -t*M*M ;
3906 int n_8 = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) -t*M*M;
3907
3908
3909 int n_9 = 2*(r) + 2*P2M * (s) +2*P2M*P2M * (t) +P2M*P2M -s*M -t*M*M;
3910 int n_10 = 2*(r) +1 + 2*P2M * (s) +2*P2M*P2M * (t) +P2M*P2M -s*M -t*M*M;
3911 int n_11 = 2*(r+1) + 2*P2M * (s) +2*P2M*P2M * (t) +P2M*P2M -s*M -t*M*M;
3912
3913 int n_12 = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t) +P2M*P2M -s*M -r -t*M*M;
3914
3915 int n_14 = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t) +P2M*P2M -s*M -(r+1) -t*M*M ;
3916
3917 int n_15 = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M -(s+1)*M -t*M*M ;
3918 int n_16 = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M -(s+1)*M -t*M*M;
3919 int n_17 = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t) +P2M*P2M -(s+1)*M -t*M*M;
3920
3921
3922 int n_18 = 2*(r) + 2*P2M * (s) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3923 int n_19 = 2*(r) +1 + 2*P2M * (s) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3924 int n_20 = 2*(r+1) + 2*P2M * (s) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3925
3926 int n_21 = 2*(r) + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3927 int n_22 = 2*(r) +1 + 2*P2M * (s) +P2M + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3928 int n_23 = 2*(r+1) + 2*P2M * (s)+P2M + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3929
3930 int n_24 = 2*(r) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3931 int n_25 = 2*(r) +1 + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3932 int n_26 = 2*(r+1) + 2*P2M * (s+1) + 2*P2M*P2M * (t+1) -(t+1)*M*M;
3933
3934
3935 // Ensuring checkerboard ordering
3936 if(((r+M*offset_x)+(s+M*offset_y)+(t+M*offset_z))%2==0){
3937 (*elementsVec)[counter] = {n_2, n_8, n_6, n_26, n_5, n_7, n_4, n_14, n_17, n_16};
3938 counter ++;
3939 (*elementsVec)[counter] = {n_2, n_6, n_0, n_18, n_4, n_3, n_1, n_10, n_12, n_9};
3940 counter ++;
3941 (*elementsVec)[counter] = {n_2, n_18, n_20, n_26, n_10, n_19, n_11, n_14, n_22, n_23};
3942 counter ++;
3943 (*elementsVec)[counter] = {n_6, n_24, n_18, n_26, n_15, n_21, n_12, n_16, n_25, n_22};
3944 counter ++;
3945 (*elementsVec)[counter] = {n_2, n_26, n_6, n_18, n_14, n_16, n_4, n_10, n_22, n_12};
3946 counter++;
3947 }
3948 else{
3949 (*elementsVec)[counter] = {n_0, n_24, n_18, n_20, n_12, n_21, n_9, n_10, n_22, n_19};
3950 counter ++;
3951 (*elementsVec)[counter] = {n_2, n_8, n_0, n_20, n_5, n_4, n_1, n_11, n_14, n_10};
3952 counter ++;
3953 (*elementsVec)[counter] = {n_8, n_6, n_0, n_24, n_7, n_3, n_4, n_16, n_15, n_12};
3954 counter ++;
3955 (*elementsVec)[counter] = {n_8, n_24, n_20, n_26, n_16, n_22, n_14, n_17, n_25, n_23};
3956 counter ++;
3957 (*elementsVec)[counter] = {n_8, n_20, n_24, n_0, n_14, n_22, n_16, n_4, n_10, n_12};
3958 counter++;
3959 }
3960
3961
3962 //std::cout << " Subwürfel zu " << " r= " << r << " s=" << s << " t= " << t << " || "<<n_0 << " " << n_1 << " " << n_2<< " " << n_3<< " " << n_4<< " " << n_5<< " " << n_6<< " " << n_7<< " " << n_8<< " " << n_9<< " " << n_10<< " " << n_11<< " " << n_12<< " " << n_14<< " " << n_15<< " " << n_16<< " " << n_17<< " " << n_18<< " " << n_19 << " " << n_20<< " " << n_21<< " " << n_22<< " " << n_23<< " " << n_24<< " " << n_25 << " " << n_26 << std::endl;
3963 }
3964 }
3965
3966 }
3967 if (verbose) {
3968 std::cout << "... done !" << std::endl;
3969 }
3970 buildElementsClass(elementsVec, elementFlag);
3971
3972}
3973template <class SC, class LO, class GO, class NO>
3974void MeshStructured<SC,LO,GO,NO>::buildSurfaces(int flagsOption, std::string FEType){
3975
3976 double tol=1.e-12;
3977 int numNodesTriangle;
3978 bool skip = false;
3979 switch (this->dim_) {
3980 case 2:
3981 break;
3982
3983 case 3:
3984 switch (flagsOption) {
3985 case 0:
3986 break;
3987 case 1:
3988 break;
3989 case 2:
3990 if(FEType == "P1")
3991 numNodesTriangle=3;
3992 else if(FEType=="P2")
3993 numNodesTriangle=6;
3994 else
3995 skip = true; //TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"For flag option and discretization no surfaces are available");
3996
3997 if(!skip){
3998 for( int T =0; T< this->elementsC_->numberElements(); T++){
3999
4000 vec_int_Type nodeList = this->elementsC_->getElement(T).getVectorNodeList();
4001
4002 vec2D_LO_Type surfaceElements_vec(4,vec_LO_Type(numNodesTriangle)); // four surfaces per element
4003
4004
4005 // Face 1 Face2 Face 3 Face 4
4006 // 2 2 * * 9 * * 3 3 * * 9 * * 2 3
4007 // * * * * * * * *
4008 // * * * * * * * *
4009 // 5 6 6 7 8 5 8 7
4010 // * * * * * * * *
4011 // * * * * * * * *
4012 // 1 * * 4 * * 0 0 1 1 * * 4 * * 0
4013 if(FEType == "P1"){
4014 surfaceElements_vec[0] = {nodeList[1],nodeList[0],nodeList[2]};
4015 surfaceElements_vec[1] = {nodeList[0],nodeList[3],nodeList[2]};
4016 surfaceElements_vec[2] = {nodeList[1],nodeList[2],nodeList[3]};
4017 surfaceElements_vec[3] = {nodeList[1],nodeList[0],nodeList[3]};
4018 }
4019 else if(FEType=="P2"){
4020 surfaceElements_vec[0] = {nodeList[1],nodeList[0],nodeList[2],nodeList[4],nodeList[6],nodeList[5]};
4021 surfaceElements_vec[1] = {nodeList[0],nodeList[3],nodeList[2],nodeList[7],nodeList[9],nodeList[6]};
4022 surfaceElements_vec[2] = {nodeList[1],nodeList[2],nodeList[3],nodeList[5],nodeList[9],nodeList[8]};
4023 surfaceElements_vec[3] = {nodeList[1],nodeList[0],nodeList[3],nodeList[4],nodeList[7],nodeList[8]};
4024 }
4025
4026 for (int i=0; i<4; i++) {
4027
4028 vec_dbl_Type p1(3),p2(3),v_E(3);
4029 p1[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(0);
4030 p1[1] =this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(1);
4031 p1[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(2);
4032
4033 p2[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(0);
4034 p2[1] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(1);
4035 p2[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(2);
4036
4037 v_E[0] = p1[1]*p2[2] - p1[2]*p2[1];
4038 v_E[1] = p1[2]*p2[0] - p1[0]*p2[2];
4039 v_E[2] = p1[0]*p2[1] - p1[1]*p2[0];
4040
4041
4042 vec_dbl_Type midpoint(3);
4043
4044 // Midpoint of triangle surface
4045 midpoint[0] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(0) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(0) )/3.;
4046 midpoint[1] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(1) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(1) )/3.;
4047 midpoint[2] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(2) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(2) )/3.;
4048
4049 int flag =10.;
4050 // x=-1 Face
4051 if (midpoint.at(0) < (coorRec[0] + tol) ) {
4052 flag = 2;
4053 if(v_E[0] > 0 )
4054 flipSurface(surfaceElements_vec[i]);
4055 }
4056
4057 p1[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(0);
4058 p1[1] =this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(1);
4059 p1[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(2);
4060
4061 p2[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(0);
4062 p2[1] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(1);
4063 p2[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(2);
4064
4065 v_E[0] = p1[1]*p2[2] - p1[2]*p2[1];
4066 v_E[1] = p1[2]*p2[0] - p1[0]*p2[2];
4067 v_E[2] = p1[0]*p2[1] - p1[1]*p2[0];
4068
4069 if(flag != 10){
4070 FiniteElement feSurface( surfaceElements_vec[i], flag);
4071 if ( !this->elementsC_->getElement(T).subElementsInitialized() )
4072 this->elementsC_->getElement(T).initializeSubElements( "P2", 2 ); // only P1 for now
4073
4074 this->elementsC_->getElement(T).addSubElement( feSurface );
4075 }
4076 }
4077 }
4078 }
4079 break;
4080 case 3:
4081 int numNodesTriangle;
4082 if(FEType == "P1")
4083 numNodesTriangle=3;
4084 else if(FEType=="P2")
4085 numNodesTriangle=6;
4086 else
4087 TEUCHOS_TEST_FOR_EXCEPTION(true,std::logic_error,"For flag option and discretization no surfaces are available");
4088
4089
4090 for( int T =0; T< this->elementsC_->numberElements(); T++){
4091
4092 vec_int_Type nodeList = this->elementsC_->getElement(T).getVectorNodeList();
4093
4094 vec2D_LO_Type surfaceElements_vec(4,vec_LO_Type(numNodesTriangle)); // four surfaces per element
4095
4096
4097 // Face 1 Face2 Face 3 Face 4
4098 // 2 2 * * 9 * * 3 3 * * 9 * * 2 3
4099 // * * * * * * * *
4100 // * * * * * * * *
4101 // 5 6 6 7 8 5 8 7
4102 // * * * * * * * *
4103 // * * * * * * * *
4104 // 1 * * 4 * * 0 0 1 1 * * 4 * * 0
4105 if(FEType == "P1"){
4106 surfaceElements_vec[0] = {nodeList[1],nodeList[0],nodeList[2]};
4107 surfaceElements_vec[1] = {nodeList[0],nodeList[3],nodeList[2]};
4108 surfaceElements_vec[2] = {nodeList[1],nodeList[2],nodeList[3]};
4109 surfaceElements_vec[3] = {nodeList[1],nodeList[0],nodeList[3]};
4110 }
4111 else if(FEType=="P2"){
4112 surfaceElements_vec[0] = {nodeList[1],nodeList[0],nodeList[2],nodeList[4],nodeList[6],nodeList[5]};
4113 surfaceElements_vec[1] = {nodeList[0],nodeList[3],nodeList[2],nodeList[7],nodeList[9],nodeList[6]};
4114 surfaceElements_vec[2] = {nodeList[1],nodeList[2],nodeList[3],nodeList[5],nodeList[9],nodeList[8]};
4115 surfaceElements_vec[3] = {nodeList[1],nodeList[0],nodeList[3],nodeList[4],nodeList[7],nodeList[8]};
4116 }
4117
4118 for (int i=0; i<4; i++) {
4119
4120 vec_dbl_Type p1(3),p2(3),v_E(3);
4121 p1[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(0);
4122 p1[1] =this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(1);
4123 p1[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(2);
4124
4125 p2[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(0);
4126 p2[1] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(1);
4127 p2[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(2);
4128
4129 v_E[0] = p1[1]*p2[2] - p1[2]*p2[1];
4130 v_E[1] = p1[2]*p2[0] - p1[0]*p2[2];
4131 v_E[2] = p1[0]*p2[1] - p1[1]*p2[0];
4132
4133
4134 vec_dbl_Type midpoint(3);
4135
4136 // Midpoint of triangle surface
4137 midpoint[0] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(0) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(0) )/3.;
4138 midpoint[1] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(1) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(1) )/3.;
4139 midpoint[2] = (this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) + this->pointsRep_->at(surfaceElements_vec[i][1]).at(2) +this->pointsRep_->at(surfaceElements_vec[i][2]).at(2) )/3.;
4140
4141 int flag =10.;
4142 // x=0 Face
4143 if (midpoint.at(0) < (coorRec[0] + tol) ) {
4144 flag = 1;
4145 if(v_E[0] > 0 )
4146 flipSurface(surfaceElements_vec[i]);
4147 }
4148 // z=0 Face
4149 if (midpoint.at(0) > (coorRec[0] + tol) &&
4150 midpoint.at(2) < (coorRec[2] + tol) ) {
4151 flag = 3;
4152 if(v_E[2] > 0 )
4153 flipSurface(surfaceElements_vec[i]);
4154 }
4155 // z=1 Face
4156 if (midpoint.at(0) > (coorRec[0] + tol) &&
4157 midpoint.at(2) > (coorRec[2] + height - tol) ) {
4158 if(v_E[2] < 0 )
4159 flipSurface(surfaceElements_vec[i]);
4160 flag = 4;
4161 }
4162 // y=0 Face
4163 if (midpoint.at(0) > (coorRec[0] + tol) &&
4164 midpoint.at(1) < (coorRec[1] + tol) ) {
4165 if(v_E[1]> 0 )
4166 flipSurface(surfaceElements_vec[i]);
4167 flag = 2;
4168 }
4169 // y=1 Face
4170 if (midpoint.at(0) > (coorRec[0] + tol) &&
4171 midpoint.at(1) > (coorRec[1] + width - tol) ) {
4172 if(v_E[1]< 0 )
4173 flipSurface(surfaceElements_vec[i]);
4174 flag = 5;
4175 }
4176 // x=1 Face
4177 if (midpoint.at(0) > (coorRec[0] + length - tol) &&
4178 midpoint.at(1) > (coorRec[1] + tol) &&
4179 midpoint.at(1) < (coorRec[1] + width - tol)&&
4180 midpoint.at(2) > (coorRec[2] + tol) &&
4181 midpoint.at(2) < (coorRec[2] + height - tol)) {
4182 if(v_E[0]< 0 )
4183 flipSurface(surfaceElements_vec[i]);
4184 flag = 6;
4185 }
4186 p1[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(0);
4187 p1[1] =this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(1);
4188 p1[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][1]).at(2);
4189
4190 p2[0] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(0) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(0);
4191 p2[1] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(1) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(1);
4192 p2[2] = this->pointsRep_->at(surfaceElements_vec[i][0]).at(2) - this->pointsRep_->at(surfaceElements_vec[i][2]).at(2);
4193
4194 v_E[0] = p1[1]*p2[2] - p1[2]*p2[1];
4195 v_E[1] = p1[2]*p2[0] - p1[0]*p2[2];
4196 v_E[2] = p1[0]*p2[1] - p1[1]*p2[0];
4197
4198 if(flag != 10){
4199 FiniteElement feSurface( surfaceElements_vec[i], flag);
4200 if ( !this->elementsC_->getElement(T).subElementsInitialized() )
4201 this->elementsC_->getElement(T).initializeSubElements( "P2", 2 ); // only P1 for now
4202
4203 this->elementsC_->getElement(T).addSubElement( feSurface );
4204 }
4205 }
4206 }
4207 break;
4208 default:
4209 break;
4210 }
4211
4212 }
4213
4214}
4215
4216// We allways want a outward normal direction
4217template <class SC, class LO, class GO, class NO>
4218void MeshStructured<SC,LO,GO,NO>::flipSurface(vec_int_Type &surfaceElements_vec){
4219
4220 LO id1,id2,id3,id4,id5,id6;
4221 id1= surfaceElements_vec[0];
4222 id2= surfaceElements_vec[1];
4223 id3= surfaceElements_vec[2];
4224 id4= surfaceElements_vec[3];
4225 id5= surfaceElements_vec[4];
4226 id6= surfaceElements_vec[5];
4227
4228 surfaceElements_vec[0] = id1;
4229 surfaceElements_vec[1] = id3;
4230 surfaceElements_vec[2] = id2;
4231 surfaceElements_vec[3] = id6;
4232 surfaceElements_vec[4] = id5;
4233 surfaceElements_vec[5] = id4;
4234}
4235
4236template <class SC, class LO, class GO, class NO>
4237void MeshStructured<SC,LO,GO,NO>::buildElementMap(){
4238
4239 Teuchos::Array<GO> elementsGlobalMapping( this->elementsC_->numberElements() );
4240 LO offset = this->comm_->getRank() * elementsGlobalMapping.size();
4241 for (int i=0; i<elementsGlobalMapping.size(); i++)
4242 elementsGlobalMapping[i] = i + offset;
4243
4244 this->elementMap_.reset(new Map<LO,GO,NO>( (GO) -1, elementsGlobalMapping(), 0, this->comm_) );
4245
4246}
4247
4248}
4249#endif
FEDD::Elements
Definition Elements.hpp:23
FEDD::FiniteElement
Definition FiniteElement.hpp:17
FEDD::Map
Definition Map_decl.hpp:30
FEDD::MeshStructured::buildMesh3DBFS
void buildMesh3DBFS(std::string FEType, int N, int M, int numProcsCoarseSolve=0)
Building general 3D backward facing step with length, width and height as defines by 'setGeometry3DBo...
Definition MeshStructured_def.hpp:2456
FEDD::MeshStructured::buildElementMap
void buildElementMap()
Building element map.
Definition MeshStructured_def.hpp:4237
FEDD::MeshStructured::build3DQ1BFS
void build3DQ1BFS(int N, int M, int numProcsCoarseSolve)
Building general 3D backward facing step with length, width and height as defines by 'setGeometry3DBo...
Definition MeshStructured_def.hpp:1627
FEDD::MeshStructured::buildMesh3D
void buildMesh3D(std::string FEType, int N, int M, int numProcsCoarseSolve=0)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox'....
Definition MeshStructured_def.hpp:575
FEDD::MeshStructured::buildP1_Disc_Q2_3DCube
void buildP1_Disc_Q2_3DCube(int N, int M, int numProcsCoarseSolve)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox' with Q2-P1 ...
Definition MeshStructured_def.hpp:964
FEDD::MeshStructured::buildMesh3D5Elements
void buildMesh3D5Elements(std::string FEType, int N, int M, int numProcsCoarseSolve)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox'....
Definition MeshStructured_def.hpp:3681
FEDD::MeshStructured::buildMesh2DBFS
void buildMesh2DBFS(std::string FEType, int N, int M, int numProcsCoarseSolve=0)
Building 2D backward facing step geometry with the lenght and height as defined per 'setGeomerty2DRec...
Definition MeshStructured_def.hpp:1999
FEDD::MeshStructured::build3DQ2_20Cube
void build3DQ2_20Cube(int N, int M, int numProcsCoarseSolve)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox' with Q2 dis...
Definition MeshStructured_def.hpp:1477
FEDD::MeshStructured::build3DQ1Cube
void build3DQ1Cube(int N, int M, int numProcsCoarseSolve)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox' with Q1 dis...
Definition MeshStructured_def.hpp:1144
FEDD::MeshStructured::buildP1_Disc_Q2_3DBFS
void buildP1_Disc_Q2_3DBFS(int N, int M, int numProcsCoarseSolve)
Building general 3D backward facing step with length, width and height as defines by 'setGeometry3DBo...
Definition MeshStructured_def.hpp:2920
FEDD::MeshStructured::build3DQ2Cube
void build3DQ2Cube(int N, int M, int numProcsCoarseSolve)
Building general 3D cuboid with length, width and height as defines by 'setGeometry3DBox' with Q2 dis...
Definition MeshStructured_def.hpp:1267
FEDD::MeshStructured::buildSurfaces
void buildSurfaces(int flagsOption, std::string FEType)
Building surfaces. This is useful for structural problems. Each surface gets another flag....
Definition MeshStructured_def.hpp:3974
FEDD::MeshStructured::buildSurfaceLinesSquare
void buildSurfaceLinesSquare()
Building suface lines aka edges. Empty.
Definition MeshStructured_def.hpp:227
FEDD::MeshStructured::setStructuredMeshFlags
void setStructuredMeshFlags(int flagsOption, std::string FEType="P1")
Setting corresponding flags to structured mesh depending on underlying problem. Rectangles/Cuboids ar...
Definition MeshStructured_def.hpp:3103
FEDD::MeshStructured::build3DQ2BFS
void build3DQ2BFS(int N, int M, int numProcsCoarseSolve)
Building general 3D backward facing step with length, width and height as defines by 'setGeometry3DBo...
Definition MeshStructured_def.hpp:1803
FEDD::MeshStructured::buildMesh2D
void buildMesh2D(std::string FEType, int N, int M, int numProcsCoarseSolve=0)
Building a general 2D rectangular mesh with the lenght and height as defined per 'setGeomerty2DRectan...
Definition MeshStructured_def.hpp:237
FEDD::Mesh
Definition Mesh_decl.hpp:25
FEDD
Adaptive Mesh Refinement.
Definition AdaptiveMeshRefinement_decl.hpp:36