# Contents of /branches/domexper/dudley/src/Mesh_tet4.c

Revision 3220 - (show annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 8 months ago) by jfenwick
File MIME type: text/plain
File size: 24897 byte(s)
```Removing references to Quadrature.? and ShapeFns

```
 1 2 /******************************************************* 3 * 4 * Copyright (c) 2003-2010 by University of Queensland 5 * Earth Systems Science Computational Center (ESSCC) 6 * http://www.uq.edu.au/esscc 7 * 8 * Primary Business: Queensland, Australia 9 * Licensed under the Open Software License version 3.0 10 * http://www.opensource.org/licenses/osl-3.0.php 11 * 12 *******************************************************/ 13 14 15 /**************************************************************/ 16 17 /* Dudley: generates rectangular meshes */ 18 19 /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */ 20 /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */ 21 /* integration scheme. */ 22 23 /**************************************************************/ 24 25 #include "TriangularMesh.h" 26 27 28 /* Be careful reading this function. The X? and NStride? are 1,2,3 but the loop vars are 0,1,2 */ 29 Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements, 30 double* Length, 31 index_t order, 32 index_t reduced_order, 33 bool_t optimize) 34 { 35 #define N_PER_E 1 36 #define DIM 3 37 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0,Nstride2=0, local_NE0, local_NE1, local_NE2, local_N0=0, local_N1=0, local_N2=0; 38 dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN; 39 index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2; 40 // Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL; 41 Dudley_Mesh* out; 42 Paso_MPIInfo *mpi_info = NULL; 43 char name[50]; 44 #ifdef Dudley_TRACE 45 double time0=Dudley_timer(); 46 #endif 47 48 const int LEFTTAG=1; /* boundary x1=0 */ 49 const int RIGHTTAG=2; /* boundary x1=1 */ 50 const int BOTTOMTAG=100; /* boundary x3=1 */ 51 const int TOPTAG=200; /* boundary x3=0 */ 52 const int FRONTTAG=10; /* boundary x2=0 */ 53 const int BACKTAG=20; /* boundary x2=1 */ 54 55 56 57 58 /* get MPI information */ 59 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD ); 60 if (! Dudley_noError()) { 61 return NULL; 62 } 63 myRank=mpi_info->rank; 64 65 /* set up the global dimensions of the mesh */ 66 67 NE0=MAX(1,numElements[0]); 68 NE1=MAX(1,numElements[1]); 69 NE2=MAX(1,numElements[2]); 70 N0=N_PER_E*NE0+1; 71 N1=N_PER_E*NE1+1; 72 N2=N_PER_E*NE2+1; 73 74 /* allocate mesh: */ 75 sprintf(name,"Triangular %d x %d x %d (x 5) mesh",N0,N1,N2); 76 out=Dudley_Mesh_alloc(name,DIM, mpi_info); 77 if (! Dudley_noError()) { 78 Paso_MPIInfo_free( mpi_info ); 79 return NULL; 80 } 81 // refElements= Dudley_ReferenceElementSet_alloc(Tet4,order,reduced_order); 82 // refFaceElements=Dudley_ReferenceElementSet_alloc(Tri3, order, reduced_order); 83 // refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order); 84 85 86 if ( Dudley_noError()) { 87 88 Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(Point1, mpi_info)); 89 Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(Tri3, mpi_info)); 90 Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(Tet4, mpi_info)); 91 92 /* work out the largest dimension */ 93 if (N2==MAX3(N0,N1,N2)) { 94 Nstride0=1; 95 Nstride1=N0; 96 Nstride2=N0*N1; 97 local_NE0=NE0; 98 e_offset0=0; 99 local_NE1=NE1; 100 e_offset1=0; 101 Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2); 102 } else if (N1==MAX3(N0,N1,N2)) { 103 Nstride0=N2; 104 Nstride1=N0*N2; 105 Nstride2=1; 106 local_NE0=NE0; 107 e_offset0=0; 108 Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1); 109 local_NE2=NE2; 110 e_offset2=0; 111 } else { 112 Nstride0=N1*N2; 113 Nstride1=1; 114 Nstride2=N1; 115 Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0); 116 local_NE1=NE1; 117 e_offset1=0; 118 local_NE2=NE2; 119 e_offset2=0; 120 } 121 offset0=e_offset0*N_PER_E; 122 offset1=e_offset1*N_PER_E; 123 offset2=e_offset2*N_PER_E; 124 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0; 125 local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0; 126 local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0; 127 128 /* get the number of surface elements */ 129 130 NFaceElements=0; 131 if (local_NE2>0) { 132 NDOF2=N2; 133 if (offset2==0) NFaceElements+=2*local_NE1*local_NE0; // each face is split 134 if (local_NE2+e_offset2 == NE2) NFaceElements+=2*local_NE1*local_NE0; 135 } else { 136 NDOF2=N2-1; 137 } 138 139 if (local_NE0>0) { 140 NDOF0=N0; 141 if (e_offset0 == 0) NFaceElements+=2*local_NE1*local_NE2; 142 if (local_NE0+e_offset0 == NE0) NFaceElements+=2*local_NE1*local_NE2; 143 } else { 144 NDOF0=N0-1; 145 } 146 147 if (local_NE1>0) { 148 NDOF1=N1; 149 if (e_offset1 == 0) NFaceElements+=2*local_NE0*local_NE2; 150 if (local_NE1+e_offset1 == NE1) NFaceElements+=2*local_NE0*local_NE2; 151 } else { 152 NDOF1=N1-1; 153 } 154 } 155 156 157 /* allocate tables: */ 158 if (Dudley_noError()) { 159 160 Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2); 161 /* we split the rectangular prism this code used to produce into 5 tetrahedrons */ 162 Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2*5); 163 /* each border face will be split in half */ 164 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements); 165 } 166 167 if (Dudley_noError()) { 168 /* create nodes */ 169 170 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2) 171 for (i2=0;i2Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0]; 179 out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1]; 180 out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2]; 181 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2; 182 out->Nodes->Tag[k]=0; 183 out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0) 184 +Nstride1*(global_i1%NDOF1) 185 +Nstride2*(global_i2%NDOF2); 186 // printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)], 187 // out->Nodes->Coordinates[INDEX2(1,k,DIM)], 188 // out->Nodes->Coordinates[INDEX2(2,k,DIM)]); 189 } 190 } 191 } 192 // printf("Now listing elements\n"); 193 /* set the elements: */ 194 195 int global_adjustment=(offset0+offset1+offset2)%2; // If we are not the only rank we may need to shift our pattern to match neighbours 196 197 NN=out->Elements->numNodes; 198 #pragma omp parallel for private(i0,i1,i2,k,node0) 199 for (i2=0;i2Elements->Id[k+j]=res+j; 210 out->Elements->Tag[k+j]=0; 211 out->Elements->Owner[k+j]=myRank; 212 } 213 214 index_t v[8]; 215 /* in non-rotated orientation the points are numbered as follows: 216 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/ 217 if ((global_adjustment+i0+i1+i2)%2==0) 218 { 219 // printf("Type0 %d, %d, %d\n",i0,i1,i2); 220 v[0]=node0; 221 v[1]=node0+Nstride0; 222 v[2]=node0+Nstride1; 223 v[3]=node0+Nstride1+Nstride0; 224 v[4]=node0+Nstride2; 225 v[5]=node0+Nstride0+Nstride2; 226 v[6]=node0+Nstride1+Nstride2; 227 v[7]=node0+Nstride2+Nstride1+Nstride0; 228 } 229 else 230 { 231 // printf("Type1 %d, %d, %d\n",i0,i1,i2); 232 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees 233 234 v[0]=node0+Nstride1; // node 0 ends up in position 2 235 v[2]=node0+Nstride1+Nstride2; // node 2 ends up in position 6 236 v[6]=node0+Nstride2; // node 6 ends up in position 4 237 v[4]=node0; // node 4 ends up in position 0 238 239 v[1]=node0+Nstride0+Nstride1; // node 1 -> pos 3 240 v[3]=node0+Nstride2+Nstride1+Nstride0; // node 3-> pos 7 241 v[7]=node0+Nstride0+Nstride2; // node 7 -> pos 5 242 v[5]=node0+Nstride0; // node 5 -> pos 1 243 } 244 245 // for (int z=0;z<8;++z) 246 // { 247 // printf("z[%d]=%d\n", z, v[z]); 248 // 249 // } 250 251 252 // index_t a=node0, b=node0+Nstride0, c=node0+Nstride1+Nstride0, d=node0+Nstride1; 253 // index_t e=node0+Nstride2, f=node0+Nstride2+Nstride0, g=node0+Nstride2+Nstride1+Nstride0, 254 // h=node0+Nstride2+Nstride1; 255 // 256 // 257 // 258 // a=0, b=1, c=3, d=2 259 // e=4, f=5, g=7, h=6 260 261 // elements nodes are numbered: centre, x, y, z 262 263 out->Elements->Nodes[INDEX2(0,k,NN)] =v[4]; 264 out->Elements->Nodes[INDEX2(1,k,NN)] =v[5]; 265 out->Elements->Nodes[INDEX2(2,k,NN)] =v[6]; 266 out->Elements->Nodes[INDEX2(3,k,NN)] =v[0]; 267 268 out->Elements->Nodes[INDEX2(0,k+1,NN)] =v[7]; 269 out->Elements->Nodes[INDEX2(1,k+1,NN)] =v[6]; 270 out->Elements->Nodes[INDEX2(2,k+1,NN)] =v[5]; 271 out->Elements->Nodes[INDEX2(3,k+1,NN)] =v[3]; 272 273 out->Elements->Nodes[INDEX2(0,k+2,NN)] =v[2]; 274 out->Elements->Nodes[INDEX2(1,k+2,NN)] =v[3]; 275 out->Elements->Nodes[INDEX2(2,k+2,NN)] =v[0]; 276 out->Elements->Nodes[INDEX2(3,k+2,NN)] =v[6]; 277 278 279 out->Elements->Nodes[INDEX2(0,k+3,NN)] =v[1]; 280 out->Elements->Nodes[INDEX2(1,k+3,NN)] =v[0]; 281 out->Elements->Nodes[INDEX2(2,k+3,NN)] =v[3]; 282 out->Elements->Nodes[INDEX2(3,k+3,NN)] =v[5]; 283 284 // I can't work out where the center is for this one 285 out->Elements->Nodes[INDEX2(0,k+4,NN)] =v[5]; 286 out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0]; 287 out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6]; 288 out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3]; 289 290 291 // for (int z=0;z<5;++z) 292 // { 293 // printf("E %d:",z); 294 // for (int q=0;q<4;++q) 295 // { 296 // index_t id=out->Elements->Nodes[INDEX2(q,z,NN)]; 297 // printf(" %d = %f, %f, %f\n", id, out->Nodes->Coordinates[INDEX2(0,id,DIM)], 298 // out->Nodes->Coordinates[INDEX2(1,id,DIM)], out->Nodes->Coordinates[INDEX2(2,id,DIM)]); 299 // } 300 // printf("\n");} 301 302 303 /* 304 for (int j=0;j<4;++j) 305 { 306 307 printf("Elt %d",j); 308 for (int m=0;m<4;++m) 309 { 310 printf(" %d",out->Elements->Nodes[INDEX2(m,k+j,NN)]); 311 } 312 printf("\n"); 313 }*/ 314 315 316 317 } 318 319 } 320 } /* end for */ 321 /* face elements */ 322 // printf("Starting face elements\n"); 323 NN=out->FaceElements->numNodes; 324 totalNECount=5*NE0*NE1*NE2; 325 faceNECount=0; 326 //printf("Bottom\n"); 327 /* these are the quadrilateral elements on boundary 1 (x3=0): */ 328 if (local_NE2>0) 329 { 330 /* ** elements on boundary 100 (x3=0): */ 331 if (e_offset2==0) 332 { 333 #pragma omp parallel for private(i0,i1,k,node0) 334 for (i1=0;i1FaceElements->Id[k]=res; 342 out->FaceElements->Tag[k]=BOTTOMTAG; 343 out->FaceElements->Owner[k]=myRank; 344 out->FaceElements->Id[k+1]=res+1; 345 out->FaceElements->Tag[k+1]=BOTTOMTAG; 346 out->FaceElements->Owner[k+1]=myRank; 347 348 index_t n0=node0; 349 index_t n1=node0+Nstride0; 350 index_t n2=node0+Nstride1; 351 index_t n3=node0+Nstride0+Nstride1; 352 353 if ((global_adjustment+i0+i1)%2==0) 354 { 355 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 356 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3; 357 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1; 358 359 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0; 360 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2; 361 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3; 362 363 } 364 else 365 { 366 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 367 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n2; 368 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n1; 369 370 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1; 371 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n2; 372 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3; 373 374 375 376 } 377 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 378 // out->FaceElements->Nodes[INDEX2(2,k,NN)]); 379 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 380 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]); 381 382 } 383 } 384 // printf("\n"); 385 faceNECount+=2*local_NE1*local_NE0; 386 } 387 totalNECount+=2*NE1*NE0; 388 389 //printf("Top\n"); 390 /* ** elements on boundary 200 (x3=1) */ 391 if (local_NE2+e_offset2 == NE2) { 392 #pragma omp parallel for private(i0,i1,k,node0) 393 for (i1=0;i1FaceElements->Id[k]=res; 401 out->FaceElements->Tag[k]=TOPTAG; 402 out->FaceElements->Owner[k]=myRank; 403 out->FaceElements->Id[k+1]=res+1; 404 out->FaceElements->Tag[k+1]=TOPTAG; 405 out->FaceElements->Owner[k+1]=myRank; 406 407 index_t n4=node0+Nstride2; 408 index_t n5=node0+Nstride0+Nstride2; 409 index_t n6=node0+Nstride1+Nstride2; 410 index_t n7=node0+Nstride1+Nstride0+Nstride2; 411 412 if ((global_adjustment+i0+i1+local_NE2-1)%2==0) 413 { 414 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4; 415 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5; 416 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6; 417 418 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n5; 419 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7; 420 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6; 421 } 422 else 423 { 424 425 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4; 426 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5; 427 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7; 428 429 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4; 430 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7; 431 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6; 432 } 433 434 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 435 // out->FaceElements->Nodes[INDEX2(2,k,NN)]); 436 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 437 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]); 438 439 } 440 } 441 /*printf("\n");*/ 442 faceNECount+=2*local_NE1*local_NE0; 443 } 444 totalNECount+=2*NE1*NE0; 445 } 446 //printf("Left\n"); 447 if (local_NE0>0) { 448 /* ** elements on boundary 001 (x1=0): */ 449 450 if (e_offset0 == 0) { 451 #pragma omp parallel for private(i1,i2,k,node0) 452 for (i2=0;i2FaceElements->Id[k]=res; 461 out->FaceElements->Tag[k]=LEFTTAG; 462 out->FaceElements->Owner[k]=myRank; 463 out->FaceElements->Id[k+1]=res+1; 464 out->FaceElements->Tag[k+1]=LEFTTAG; 465 out->FaceElements->Owner[k+1]=myRank; 466 467 index_t n0=node0; 468 index_t n2=node0+Nstride1; 469 index_t n4=node0+Nstride2; 470 index_t n6=node0+Nstride1+Nstride2; 471 472 if ((global_adjustment+0+i1+i2)%2==0) 473 { 474 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 475 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4; 476 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6; 477 478 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0; 479 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6; 480 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2; 481 } 482 else 483 { 484 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees 485 486 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 487 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4; 488 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2; 489 490 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4; 491 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6; 492 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2; 493 } 494 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 495 // out->FaceElements->Nodes[INDEX2(2,k,NN)]); 496 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 497 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]); 498 } 499 } 500 faceNECount+=2*local_NE1*local_NE2; 501 } 502 /*printf("\n");*/ 503 //printf("Right\n"); 504 totalNECount+=2*NE1*NE2; 505 /* ** elements on boundary 002 (x1=1): */ 506 if (local_NE0+e_offset0 == NE0) { 507 #pragma omp parallel for private(i1,i2,k,node0) 508 for (i2=0;i2FaceElements->Id[k]=res; 516 out->FaceElements->Tag[k]=RIGHTTAG; 517 out->FaceElements->Owner[k]=myRank; 518 out->FaceElements->Id[k+1]=res+1; 519 out->FaceElements->Tag[k+1]=RIGHTTAG; 520 out->FaceElements->Owner[k+1]=myRank; 521 522 index_t n1=node0+Nstride0; 523 index_t n3=node0+Nstride0+Nstride1; 524 index_t n5=node0+Nstride0+Nstride2; 525 index_t n7=node0+Nstride0+Nstride1+Nstride2; 526 if ((global_adjustment+local_NE0-1+i1+i2)%2==0) 527 { 528 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1; 529 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3; 530 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5; 531 532 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n3; 533 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7; 534 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n5; 535 } 536 else 537 { 538 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees 539 540 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1; 541 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n7; 542 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5; 543 544 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1; 545 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3; 546 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n7; 547 } 548 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 549 // out->FaceElements->Nodes[INDEX2(2,k,NN)]); 550 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 551 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]); 552 553 } 554 } 555 faceNECount+=2*local_NE1*local_NE2; 556 } 557 totalNECount+=2*NE1*NE2; 558 } 559 //printf("\n"); 560 //printf("Front\n"); 561 if (local_NE1>0) { 562 /* ** elements on boundary 010 (x2=0): */ 563 if (e_offset1 == 0) { 564 #pragma omp parallel for private(i0,i2,k,node0) 565 for (i2=0;i2FaceElements->Id[k]=res; 572 out->FaceElements->Tag[k]=FRONTTAG; 573 out->FaceElements->Owner[k]=myRank; 574 out->FaceElements->Id[k+1]=res+1; 575 out->FaceElements->Tag[k+1]=FRONTTAG; 576 out->FaceElements->Owner[k+1]=myRank; 577 578 index_t n0=node0; 579 index_t n1=node0+Nstride0; 580 index_t n4=node0+Nstride2; 581 index_t n5=node0+Nstride0+Nstride2; 582 583 if ((global_adjustment+i0+0+i2)%2==0) 584 { 585 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 586 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1; 587 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5; 588 589 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0; 590 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5; 591 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4; 592 593 594 } 595 else 596 { 597 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees 598 599 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0; 600 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1; 601 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n4; 602 603 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1; 604 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5; 605 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4; 606 607 } 608 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 609 out->FaceElements->Nodes[INDEX2(2,k,NN)]); 610 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 611 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/ 612 } 613 } 614 faceNECount+=2*local_NE0*local_NE2; 615 } 616 // printf("\n"); 617 totalNECount+=2*NE0*NE2; 618 //printf("Back\n"); 619 /* ** elements on boundary 020 (x2=1): */ 620 if (local_NE1+e_offset1 == NE1) { 621 #pragma omp parallel for private(i0,i2,k,node0) 622 for (i2=0;i2FaceElements->Id[k]=res; 629 out->FaceElements->Tag[k]=BACKTAG; 630 out->FaceElements->Owner[k]=myRank; 631 out->FaceElements->Id[k+1]=res+1; 632 out->FaceElements->Tag[k+1]=BACKTAG; 633 out->FaceElements->Owner[k+1]=myRank; 634 635 index_t n2=node0+Nstride1; 636 index_t n6=node0+Nstride1+Nstride2; 637 index_t n7=node0+Nstride0+Nstride1+Nstride2; 638 index_t n3=node0+Nstride0+Nstride1; 639 640 if ((global_adjustment+i0+local_NE1-1+i2)%2==0) 641 { 642 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2; 643 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6; 644 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n3; 645 646 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n6; 647 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7; 648 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3; 649 650 651 } 652 else 653 { 654 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees 655 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2; 656 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6; 657 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7; 658 659 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n2; 660 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7; 661 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3; 662 663 } 664 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)], 665 out->FaceElements->Nodes[INDEX2(2,k,NN)]); 666 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)], 667 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/ 668 } 669 } 670 faceNECount+=2*local_NE0*local_NE2; 671 } 672 totalNECount+=2*NE0*NE2; 673 } 674 } 675 // printf("\n"); 676 if (Dudley_noError()) { 677 /* add tag names */ 678 Dudley_Mesh_addTagMap(out,"top", TOPTAG); 679 Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG); 680 Dudley_Mesh_addTagMap(out,"left", LEFTTAG); 681 Dudley_Mesh_addTagMap(out,"right", RIGHTTAG); 682 Dudley_Mesh_addTagMap(out,"front", FRONTTAG); 683 Dudley_Mesh_addTagMap(out,"back", BACKTAG); 684 } 685 /* prepare mesh for further calculatuions:*/ 686 if (Dudley_noError()) { 687 Dudley_Mesh_resolveNodeIds(out); 688 } 689 if (Dudley_noError()) { 690 Dudley_Mesh_prepare(out, optimize); 691 } 692 693 if (!Dudley_noError()) { 694 Dudley_Mesh_free(out); 695 } 696 /* free up memory */ 697 // Dudley_ReferenceElementSet_dealloc(refPoints); 698 // Dudley_ReferenceElementSet_dealloc(refFaceElements); 699 // Dudley_ReferenceElementSet_dealloc(refElements); 700 Paso_MPIInfo_free( mpi_info ); 701 702 return out; 703 } 704 705

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo

 ViewVC Help Powered by ViewVC 1.1.26