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