/[escript]/branches/domexper/dudley/src/Mesh_tet4.c
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo /branches/scons_revamp_from_3210/dudley/src/Mesh_tet4.c:3212-3243

  ViewVC Help
Powered by ViewVC 1.1.26