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