# Contents of /branches/trilinos_from_5897/dudley/src/Mesh_tet4.cpp

Revision 6009 - (show annotations)
Wed Mar 2 04:13:26 2016 UTC (3 years, 1 month ago) by caltinay
File size: 25614 byte(s)
```Much needed sync with trunk...

```
 1 2 /***************************************************************************** 3 * 4 * Copyright (c) 2003-2016 by The University of Queensland 5 6 * 7 * Primary Business: Queensland, Australia 8 * Licensed under the Open Software License version 3.0 9 10 * 11 * Development until 2012 by Earth Systems Science Computational Center (ESSCC) 12 * Development 2012-2013 by School of Earth Sciences 13 * Development from 2014 by Centre for Geoscience Computing (GeoComp) 14 * 15 *****************************************************************************/ 16 17 /****************************************************************************/ 18 19 /* Dudley: generates rectangular meshes */ 20 21 /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */ 22 /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */ 23 /* integration scheme. */ 24 25 /****************************************************************************/ 26 27 #include "TriangularMesh.h" 28 29 #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_)) 30 31 namespace dudley { 32 33 // Be careful reading this function. The X? and NStride? are 1,2,3 34 // but the loop vars are 0,1,2 35 Dudley_Mesh *Dudley_TriangularMesh_Tet4(dim_t * numElements, 36 double *Length, index_t order, index_t reduced_order, bool optimize, escript::JMPI& mpi_info) 37 { 38 #define N_PER_E 1 39 #define DIM 3 40 dim_t N0, N1, N2, NE0, NE1, NE2, i0, i1, i2, k, Nstride0 = 0, Nstride1 = 0, Nstride2 = 41 0, local_NE0, local_NE1, local_NE2, local_N0 = 0, local_N1 = 0, local_N2 = 0; 42 dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NDOF2 = 0, NFaceElements = 0, NN; 43 index_t node0, myRank, e_offset2, e_offset1, e_offset0 = 0, offset1 = 0, offset2 = 0, offset0 = 44 0, global_i0, global_i1, global_i2; 45 Dudley_Mesh *out; 46 char name[50]; 47 #ifdef Dudley_TRACE 48 double time0 = Dudley_timer(); 49 #endif 50 51 const int LEFTTAG = 1; /* boundary x1=0 */ 52 const int RIGHTTAG = 2; /* boundary x1=1 */ 53 const int BOTTOMTAG = 100; /* boundary x3=1 */ 54 const int TOPTAG = 200; /* boundary x3=0 */ 55 const int FRONTTAG = 10; /* boundary x2=0 */ 56 const int BACKTAG = 20; /* boundary x2=1 */ 57 58 /* get MPI information */ 59 myRank = mpi_info->rank; 60 61 /* set up the global dimensions of the mesh */ 62 NE0 = std::max(1, numElements[0]); 63 NE1 = std::max(1, numElements[1]); 64 NE2 = std::max(1, numElements[2]); 65 N0 = N_PER_E * NE0 + 1; 66 N1 = N_PER_E * NE1 + 1; 67 N2 = N_PER_E * NE2 + 1; 68 69 /* allocate mesh: */ 70 sprintf(name, "Triangular %d x %d x %d (x 5) mesh", N0, N1, N2); 71 out = Dudley_Mesh_alloc(name, DIM, mpi_info); 72 73 Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info)); 74 Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info)); 75 Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tet4, mpi_info)); 76 77 /* work out the largest dimension */ 78 if (N2 == MAX3(N0, N1, N2)) { 79 Nstride0 = 1; 80 Nstride1 = N0; 81 Nstride2 = N0 * N1; 82 local_NE0 = NE0; 83 e_offset0 = 0; 84 local_NE1 = NE1; 85 e_offset1 = 0; 86 mpi_info->split(NE2, &local_NE2, &e_offset2); 87 } else if (N1 == MAX3(N0, N1, N2)) { 88 Nstride0 = N2; 89 Nstride1 = N0 * N2; 90 Nstride2 = 1; 91 local_NE0 = NE0; 92 e_offset0 = 0; 93 mpi_info->split(NE1, &local_NE1, &e_offset1); 94 local_NE2 = NE2; 95 e_offset2 = 0; 96 } else { 97 Nstride0 = N1 * N2; 98 Nstride1 = 1; 99 Nstride2 = N1; 100 mpi_info->split(NE0, &local_NE0, &e_offset0); 101 local_NE1 = NE1; 102 e_offset1 = 0; 103 local_NE2 = NE2; 104 e_offset2 = 0; 105 } 106 offset0 = e_offset0 * N_PER_E; 107 offset1 = e_offset1 * N_PER_E; 108 offset2 = e_offset2 * N_PER_E; 109 local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0; 110 local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0; 111 local_N2 = local_NE2 > 0 ? local_NE2 * N_PER_E + 1 : 0; 112 113 /* get the number of surface elements */ 114 115 NFaceElements = 0; 116 if (local_NE2 > 0) { 117 NDOF2 = N2; 118 if (offset2 == 0) 119 NFaceElements += 2 * local_NE1 * local_NE0; /* each face is split */ 120 if (local_NE2 + e_offset2 == NE2) 121 NFaceElements += 2 * local_NE1 * local_NE0; 122 } else { 123 NDOF2 = N2 - 1; 124 } 125 126 if (local_NE0 > 0) { 127 NDOF0 = N0; 128 if (e_offset0 == 0) 129 NFaceElements += 2 * local_NE1 * local_NE2; 130 if (local_NE0 + e_offset0 == NE0) 131 NFaceElements += 2 * local_NE1 * local_NE2; 132 } else { 133 NDOF0 = N0 - 1; 134 } 135 136 if (local_NE1 > 0) { 137 NDOF1 = N1; 138 if (e_offset1 == 0) 139 NFaceElements += 2 * local_NE0 * local_NE2; 140 if (local_NE1 + e_offset1 == NE1) 141 NFaceElements += 2 * local_NE0 * local_NE2; 142 } else { 143 NDOF1 = N1 - 1; 144 } 145 146 /* allocate tables: */ 147 Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1 * local_N2); 148 /* we split the rectangular prism this code used to produce into 5 tetrahedrons */ 149 Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * local_NE2 * 5); 150 /* each border face will be split in half */ 151 Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements); 152 153 int global_adjustment; 154 155 /* create nodes */ 156 157 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2) 158 for (i2 = 0; i2 < local_N2; i2++) { 159 for (i1 = 0; i1 < local_N1; i1++) { 160 for (i0 = 0; i0 < local_N0; i0++) { 161 k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2; 162 global_i0 = i0 + offset0; 163 global_i1 = i1 + offset1; 164 global_i2 = i2 + offset2; 165 out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0]; 166 out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1]; 167 out->Nodes->Coordinates[INDEX2(2, k, DIM)] = (double)global_i2 / (double)(N2 - 1) * Length[2]; 168 out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2; 169 out->Nodes->Tag[k] = 0; 170 out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0) 171 + Nstride1 * (global_i1 % NDOF1) + Nstride2 * (global_i2 % NDOF2); 172 } 173 } 174 } 175 /* set the elements: */ 176 177 global_adjustment = (offset0 + offset1 + offset2) % 2; /* If we are not the only rank we may need to shift our pattern to match neighbours */ 178 179 NN = out->Elements->numNodes; 180 #pragma omp parallel for private(i0,i1,i2,k,node0) 181 for (i2 = 0; i2 < local_NE2; i2++) { 182 for (i1 = 0; i1 < local_NE1; i1++) { 183 for (i0 = 0; i0 < local_NE0; i0++) { 184 index_t res; 185 index_t v[8]; 186 int j; 187 k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2); 188 node0 = 189 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) + 190 Nstride2 * N_PER_E * (i2 + e_offset2); 191 192 res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2)); 193 for (j = 0; j < 5; ++j) { 194 out->Elements->Id[k + j] = res + j; 195 out->Elements->Tag[k + j] = 0; 196 out->Elements->Owner[k + j] = myRank; 197 } 198 199 /* in non-rotated orientation the points are numbered as follows: 200 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/ 201 if ((global_adjustment + i0 + i1 + i2) % 2 == 0) { 202 v[0] = node0; 203 v[1] = node0 + Nstride0; 204 v[2] = node0 + Nstride1; 205 v[3] = node0 + Nstride1 + Nstride0; 206 v[4] = node0 + Nstride2; 207 v[5] = node0 + Nstride0 + Nstride2; 208 v[6] = node0 + Nstride1 + Nstride2; 209 v[7] = node0 + Nstride2 + Nstride1 + Nstride0; 210 } else { 211 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */ 212 213 v[0] = node0 + Nstride1; /* node 0 ends up in position 2 */ 214 v[2] = node0 + Nstride1 + Nstride2; /* node 2 ends up in position 6 */ 215 v[6] = node0 + Nstride2; /* node 6 ends up in position 4 */ 216 v[4] = node0; /* node 4 ends up in position 0 */ 217 218 v[1] = node0 + Nstride0 + Nstride1; /* node 1 -> pos 3 */ 219 v[3] = node0 + Nstride2 + Nstride1 + Nstride0; /* node 3-> pos 7 */ 220 v[7] = node0 + Nstride0 + Nstride2; /* node 7 -> pos 5 */ 221 v[5] = node0 + Nstride0; /* node 5 -> pos 1 */ 222 } 223 224 /* elements nodes are numbered: centre, x, y, z */ 225 226 out->Elements->Nodes[INDEX2(0, k, NN)] = v[4]; 227 out->Elements->Nodes[INDEX2(1, k, NN)] = v[5]; 228 out->Elements->Nodes[INDEX2(2, k, NN)] = v[6]; 229 out->Elements->Nodes[INDEX2(3, k, NN)] = v[0]; 230 231 out->Elements->Nodes[INDEX2(0, k + 1, NN)] = v[7]; 232 out->Elements->Nodes[INDEX2(1, k + 1, NN)] = v[6]; 233 out->Elements->Nodes[INDEX2(2, k + 1, NN)] = v[5]; 234 out->Elements->Nodes[INDEX2(3, k + 1, NN)] = v[3]; 235 236 out->Elements->Nodes[INDEX2(0, k + 2, NN)] = v[2]; 237 out->Elements->Nodes[INDEX2(1, k + 2, NN)] = v[3]; 238 out->Elements->Nodes[INDEX2(2, k + 2, NN)] = v[0]; 239 out->Elements->Nodes[INDEX2(3, k + 2, NN)] = v[6]; 240 241 out->Elements->Nodes[INDEX2(0, k + 3, NN)] = v[1]; 242 out->Elements->Nodes[INDEX2(1, k + 3, NN)] = v[0]; 243 out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3]; 244 out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5]; 245 246 /* I can't work out where the center is for this one */ 247 out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5]; 248 out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0]; 249 out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6]; 250 out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3]; 251 } 252 } 253 } /* end for */ 254 /* face elements */ 255 NN = out->FaceElements->numNodes; 256 totalNECount = 5 * NE0 * NE1 * NE2; 257 faceNECount = 0; 258 /* these are the quadrilateral elements on boundary 1 (x3=0): */ 259 if (local_NE2 > 0) { 260 /* ** elements on boundary 100 (x3=0): */ 261 if (e_offset2 == 0) { 262 #pragma omp parallel for private(i0,i1,k,node0) 263 for (i1 = 0; i1 < local_NE1; i1++) { 264 for (i0 = 0; i0 < local_NE0; i0++) { 265 index_t res, n0, n1, n2, n3; 266 k = 2 * (i0 + local_NE0 * i1) + faceNECount; 267 node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1); 268 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount; 269 out->FaceElements->Id[k] = res; 270 out->FaceElements->Tag[k] = BOTTOMTAG; 271 out->FaceElements->Owner[k] = myRank; 272 out->FaceElements->Id[k + 1] = res + 1; 273 out->FaceElements->Tag[k + 1] = BOTTOMTAG; 274 out->FaceElements->Owner[k + 1] = myRank; 275 276 n0 = node0; 277 n1 = node0 + Nstride0; 278 n2 = node0 + Nstride1; 279 n3 = node0 + Nstride0 + Nstride1; 280 281 if ((global_adjustment + i0 + i1) % 2 == 0) { 282 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 283 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3; 284 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1; 285 286 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0; 287 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2; 288 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3; 289 290 } else { 291 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 292 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n2; 293 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1; 294 295 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1; 296 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2; 297 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3; 298 299 } 300 } 301 } 302 faceNECount += 2 * local_NE1 * local_NE0; 303 } 304 totalNECount += 2 * NE1 * NE0; 305 /* ** elements on boundary 200 (x3=1) - Top */ 306 if (local_NE2 + e_offset2 == NE2) { 307 #pragma omp parallel for private(i0,i1,k,node0) 308 for (i1 = 0; i1 < local_NE1; i1++) { 309 for (i0 = 0; i0 < local_NE0; i0++) { 310 311 index_t res, n4, n5, n6, n7; 312 k = 2 * (i0 + local_NE0 * i1) + faceNECount; 313 node0 = 314 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) + 315 Nstride2 * N_PER_E * (NE2 - 1); 316 317 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount; 318 out->FaceElements->Id[k] = res; 319 out->FaceElements->Tag[k] = TOPTAG; 320 out->FaceElements->Owner[k] = myRank; 321 out->FaceElements->Id[k + 1] = res + 1; 322 out->FaceElements->Tag[k + 1] = TOPTAG; 323 out->FaceElements->Owner[k + 1] = myRank; 324 325 n4 = node0 + Nstride2; 326 n5 = node0 + Nstride0 + Nstride2; 327 n6 = node0 + Nstride1 + Nstride2; 328 n7 = node0 + Nstride1 + Nstride0 + Nstride2; 329 330 if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0) { 331 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4; 332 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5; 333 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6; 334 335 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n5; 336 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7; 337 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6; 338 } else { 339 340 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4; 341 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5; 342 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7; 343 344 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4; 345 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7; 346 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6; 347 } 348 } 349 } 350 faceNECount += 2 * local_NE1 * local_NE0; 351 } 352 totalNECount += 2 * NE1 * NE0; 353 } 354 if (local_NE0 > 0) { 355 /* ** elements on boundary 001 (x1=0): - Left */ 356 357 if (e_offset0 == 0) { 358 #pragma omp parallel for private(i1,i2,k,node0) 359 for (i2 = 0; i2 < local_NE2; i2++) { 360 for (i1 = 0; i1 < local_NE1; i1++) { 361 362 index_t res, n0, n2, n4, n6; 363 k = 2 * (i1 + local_NE1 * i2) + faceNECount; 364 node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2); 365 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount; 366 out->FaceElements->Id[k] = res; 367 out->FaceElements->Tag[k] = LEFTTAG; 368 out->FaceElements->Owner[k] = myRank; 369 out->FaceElements->Id[k + 1] = res + 1; 370 out->FaceElements->Tag[k + 1] = LEFTTAG; 371 out->FaceElements->Owner[k + 1] = myRank; 372 373 n0 = node0; 374 n2 = node0 + Nstride1; 375 n4 = node0 + Nstride2; 376 n6 = node0 + Nstride1 + Nstride2; 377 378 if ((global_adjustment + 0 + i1 + i2) % 2 == 0) { 379 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 380 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4; 381 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6; 382 383 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0; 384 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6; 385 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2; 386 } else { 387 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */ 388 389 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 390 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4; 391 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2; 392 393 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4; 394 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6; 395 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2; 396 } 397 } 398 } 399 faceNECount += 2 * local_NE1 * local_NE2; 400 } 401 totalNECount += 2 * NE1 * NE2; 402 /* ** elements on boundary 002 (x1=1): - Right */ 403 if (local_NE0 + e_offset0 == NE0) { 404 #pragma omp parallel for private(i1,i2,k,node0) 405 for (i2 = 0; i2 < local_NE2; i2++) { 406 for (i1 = 0; i1 < local_NE1; i1++) { 407 index_t res, n1, n3, n5, n7; 408 k = 2 * (i1 + local_NE1 * i2) + faceNECount; 409 410 node0 = 411 Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) + 412 Nstride2 * N_PER_E * (i2 + e_offset2); 413 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount; 414 out->FaceElements->Id[k] = res; 415 out->FaceElements->Tag[k] = RIGHTTAG; 416 out->FaceElements->Owner[k] = myRank; 417 out->FaceElements->Id[k + 1] = res + 1; 418 out->FaceElements->Tag[k + 1] = RIGHTTAG; 419 out->FaceElements->Owner[k + 1] = myRank; 420 421 n1 = node0 + Nstride0; 422 n3 = node0 + Nstride0 + Nstride1; 423 n5 = node0 + Nstride0 + Nstride2; 424 n7 = node0 + Nstride0 + Nstride1 + Nstride2; 425 if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0) { 426 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1; 427 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3; 428 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5; 429 430 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n3; 431 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7; 432 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5; 433 } else { 434 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */ 435 436 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1; 437 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7; 438 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5; 439 440 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1; 441 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n3; 442 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n7; 443 } 444 } 445 } 446 faceNECount += 2 * local_NE1 * local_NE2; 447 } 448 totalNECount += 2 * NE1 * NE2; 449 } 450 if (local_NE1 > 0) { 451 /* ** elements on boundary 010 (x2=0): -Front */ 452 if (e_offset1 == 0) { 453 #pragma omp parallel for private(i0,i2,k,node0) 454 for (i2 = 0; i2 < local_NE2; i2++) { 455 for (i0 = 0; i0 < local_NE0; i0++) { 456 index_t res, n0, n1, n4, n5; 457 k = 2 * (i0 + local_NE0 * i2) + faceNECount; 458 node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2); 459 res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount; 460 out->FaceElements->Id[k] = res; 461 out->FaceElements->Tag[k] = FRONTTAG; 462 out->FaceElements->Owner[k] = myRank; 463 out->FaceElements->Id[k + 1] = res + 1; 464 out->FaceElements->Tag[k + 1] = FRONTTAG; 465 out->FaceElements->Owner[k + 1] = myRank; 466 467 n0 = node0; 468 n1 = node0 + Nstride0; 469 n4 = node0 + Nstride2; 470 n5 = node0 + Nstride0 + Nstride2; 471 472 if ((global_adjustment + i0 + 0 + i2) % 2 == 0) { 473 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 474 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1; 475 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5; 476 477 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0; 478 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5; 479 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4; 480 481 } else { 482 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */ 483 484 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0; 485 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1; 486 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4; 487 488 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1; 489 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5; 490 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4; 491 492 } 493 } 494 } 495 faceNECount += 2 * local_NE0 * local_NE2; 496 } 497 totalNECount += 2 * NE0 * NE2; 498 /* ** elements on boundary 020 (x2=1): - Back */ 499 if (local_NE1 + e_offset1 == NE1) { 500 #pragma omp parallel for private(i0,i2,k,node0) 501 for (i2 = 0; i2 < local_NE2; i2++) { 502 for (i0 = 0; i0 < local_NE0; i0++) { 503 index_t res, n2, n6, n7, n3; 504 k = 2 * (i0 + local_NE0 * i2) + faceNECount; 505 node0 = 506 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1) + 507 Nstride2 * N_PER_E * (i2 + e_offset2); 508 res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount; 509 out->FaceElements->Id[k] = res; 510 out->FaceElements->Tag[k] = BACKTAG; 511 out->FaceElements->Owner[k] = myRank; 512 out->FaceElements->Id[k + 1] = res + 1; 513 out->FaceElements->Tag[k + 1] = BACKTAG; 514 out->FaceElements->Owner[k + 1] = myRank; 515 516 n2 = node0 + Nstride1; 517 n6 = node0 + Nstride1 + Nstride2; 518 n7 = node0 + Nstride0 + Nstride1 + Nstride2; 519 n3 = node0 + Nstride0 + Nstride1; 520 521 if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0) { 522 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2; 523 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6; 524 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n3; 525 526 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n6; 527 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7; 528 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3; 529 530 } else { 531 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */ 532 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2; 533 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6; 534 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7; 535 536 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2; 537 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7; 538 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3; 539 540 } 541 } 542 } 543 faceNECount += 2 * local_NE0 * local_NE2; 544 } 545 totalNECount += 2 * NE0 * NE2; 546 } 547 /* add tag names */ 548 Dudley_Mesh_addTagMap(out, "top", TOPTAG); 549 Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG); 550 Dudley_Mesh_addTagMap(out, "left", LEFTTAG); 551 Dudley_Mesh_addTagMap(out, "right", RIGHTTAG); 552 Dudley_Mesh_addTagMap(out, "front", FRONTTAG); 553 Dudley_Mesh_addTagMap(out, "back", BACKTAG); 554 // prepare mesh for further calculations 555 Dudley_Mesh_resolveNodeIds(out); 556 Dudley_Mesh_prepare(out, optimize); 557 558 return out; 559 } 560 561 } // namespace dudley 562

## Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/4.0fordebian/dudley/src/Mesh_tet4.cpp:5567-5588 /branches/pasowrap/dudley/src/Mesh_tet4.c:3661-3674 /branches/py3_attempt2/dudley/src/Mesh_tet4.c:3871-3891 /branches/ripleygmg_from_3668/dudley/src/Mesh_tet4.c:3669-3791 /branches/scons_revamp_from_3210/dudley/src/Mesh_tet4.c:3212-3243 /branches/symbolic_from_3470/dudley/src/Mesh_tet4.c:3471-3974 /branches/symbolic_from_3470/ripley/test/python/dudley/src/Mesh_tet4.c:3517-3974 /release/4.0/dudley/src/Mesh_tet4.cpp:5380-5406 /trunk/dudley/src/Mesh_tet4.cpp:4257-4344,5898-6007 /trunk/ripley/test/python/dudley/src/Mesh_tet4.c:3480-3515

 ViewVC Help Powered by ViewVC 1.1.26