/[escript]/trunk/dudley/src/Mesh_tet4.c
ViewVC logotype

Contents of /trunk/dudley/src/Mesh_tet4.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3259 - (show annotations)
Mon Oct 11 01:48:14 2010 UTC (9 years ago) by jfenwick
File MIME type: text/plain
File size: 21169 byte(s)
Merging dudley and scons updates from branches

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 int global_adjustment;
181
182 /* create nodes */
183
184 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
185 for (i2 = 0; i2 < local_N2; i2++)
186 {
187 for (i1 = 0; i1 < local_N1; i1++)
188 {
189 for (i0 = 0; i0 < local_N0; i0++)
190 {
191 k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2;
192 global_i0 = i0 + offset0;
193 global_i1 = i1 + offset1;
194 global_i2 = i2 + offset2;
195 out->Nodes->Coordinates[INDEX2(0, k, DIM)] = DBLE(global_i0) / DBLE(N0 - 1) * Length[0];
196 out->Nodes->Coordinates[INDEX2(1, k, DIM)] = DBLE(global_i1) / DBLE(N1 - 1) * Length[1];
197 out->Nodes->Coordinates[INDEX2(2, k, DIM)] = DBLE(global_i2) / DBLE(N2 - 1) * Length[2];
198 out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2;
199 out->Nodes->Tag[k] = 0;
200 out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)
201 + Nstride1 * (global_i1 % NDOF1) + Nstride2 * (global_i2 % NDOF2);
202 }
203 }
204 }
205 /* set the elements: */
206
207 global_adjustment = (offset0 + offset1 + offset2) % 2; /* If we are not the only rank we may need to shift our pattern to match neighbours */
208
209 NN = out->Elements->numNodes;
210 #pragma omp parallel for private(i0,i1,i2,k,node0)
211 for (i2 = 0; i2 < local_NE2; i2++)
212 {
213 for (i1 = 0; i1 < local_NE1; i1++)
214 {
215 for (i0 = 0; i0 < local_NE0; i0++)
216 {
217 index_t res;
218 index_t v[8];
219 int j;
220 k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2);
221 node0 =
222 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
223 Nstride2 * N_PER_E * (i2 + e_offset2);
224
225 res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2));
226 for (j = 0; j < 5; ++j)
227 {
228 out->Elements->Id[k + j] = res + j;
229 out->Elements->Tag[k + j] = 0;
230 out->Elements->Owner[k + j] = myRank;
231 }
232
233 /* in non-rotated orientation the points are numbered as follows:
234 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
235 if ((global_adjustment + i0 + i1 + i2) % 2 == 0)
236 {
237 v[0] = node0;
238 v[1] = node0 + Nstride0;
239 v[2] = node0 + Nstride1;
240 v[3] = node0 + Nstride1 + Nstride0;
241 v[4] = node0 + Nstride2;
242 v[5] = node0 + Nstride0 + Nstride2;
243 v[6] = node0 + Nstride1 + Nstride2;
244 v[7] = node0 + Nstride2 + Nstride1 + Nstride0;
245 }
246 else
247 {
248 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
249
250 v[0] = node0 + Nstride1; /* node 0 ends up in position 2 */
251 v[2] = node0 + Nstride1 + Nstride2; /* node 2 ends up in position 6 */
252 v[6] = node0 + Nstride2; /* node 6 ends up in position 4 */
253 v[4] = node0; /* node 4 ends up in position 0 */
254
255 v[1] = node0 + Nstride0 + Nstride1; /* node 1 -> pos 3 */
256 v[3] = node0 + Nstride2 + Nstride1 + Nstride0; /* node 3-> pos 7 */
257 v[7] = node0 + Nstride0 + Nstride2; /* node 7 -> pos 5 */
258 v[5] = node0 + Nstride0; /* node 5 -> pos 1 */
259 }
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 out->Elements->Nodes[INDEX2(0, k + 3, NN)] = v[1];
279 out->Elements->Nodes[INDEX2(1, k + 3, NN)] = v[0];
280 out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3];
281 out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5];
282
283 /* I can't work out where the center is for this one */
284 out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5];
285 out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0];
286 out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6];
287 out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3];
288 }
289
290 }
291 } /* end for */
292 /* face elements */
293 NN = out->FaceElements->numNodes;
294 totalNECount = 5 * NE0 * NE1 * NE2;
295 faceNECount = 0;
296 /* these are the quadrilateral elements on boundary 1 (x3=0): */
297 if (local_NE2 > 0)
298 {
299 /* ** elements on boundary 100 (x3=0): */
300 if (e_offset2 == 0)
301 {
302 #pragma omp parallel for private(i0,i1,k,node0)
303 for (i1 = 0; i1 < local_NE1; i1++)
304 {
305 for (i0 = 0; i0 < local_NE0; i0++)
306 {
307 index_t res, n0, n1, n2, n3;
308 k = 2 * (i0 + local_NE0 * i1) + faceNECount;
309 node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);
310 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
311 out->FaceElements->Id[k] = res;
312 out->FaceElements->Tag[k] = BOTTOMTAG;
313 out->FaceElements->Owner[k] = myRank;
314 out->FaceElements->Id[k + 1] = res + 1;
315 out->FaceElements->Tag[k + 1] = BOTTOMTAG;
316 out->FaceElements->Owner[k + 1] = myRank;
317
318 n0 = node0;
319 n1 = node0 + Nstride0;
320 n2 = node0 + Nstride1;
321 n3 = node0 + Nstride0 + Nstride1;
322
323 if ((global_adjustment + i0 + i1) % 2 == 0)
324 {
325 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
326 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;
327 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;
328
329 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
330 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;
331 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
332
333 }
334 else
335 {
336 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
337 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n2;
338 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;
339
340 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
341 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;
342 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
343
344 }
345 }
346 }
347 faceNECount += 2 * local_NE1 * local_NE0;
348 }
349 totalNECount += 2 * NE1 * NE0;
350 /* ** elements on boundary 200 (x3=1) - Top */
351 if (local_NE2 + e_offset2 == NE2)
352 {
353 #pragma omp parallel for private(i0,i1,k,node0)
354 for (i1 = 0; i1 < local_NE1; i1++)
355 {
356 for (i0 = 0; i0 < local_NE0; i0++)
357 {
358
359 index_t res, n4, n5, n6, n7;
360 k = 2 * (i0 + local_NE0 * i1) + faceNECount;
361 node0 =
362 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
363 Nstride2 * N_PER_E * (NE2 - 1);
364
365 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
366 out->FaceElements->Id[k] = res;
367 out->FaceElements->Tag[k] = TOPTAG;
368 out->FaceElements->Owner[k] = myRank;
369 out->FaceElements->Id[k + 1] = res + 1;
370 out->FaceElements->Tag[k + 1] = TOPTAG;
371 out->FaceElements->Owner[k + 1] = myRank;
372
373 n4 = node0 + Nstride2;
374 n5 = node0 + Nstride0 + Nstride2;
375 n6 = node0 + Nstride1 + Nstride2;
376 n7 = node0 + Nstride1 + Nstride0 + Nstride2;
377
378 if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0)
379 {
380 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;
381 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;
382 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;
383
384 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n5;
385 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
386 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;
387 }
388 else
389 {
390
391 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;
392 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;
393 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;
394
395 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;
396 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
397 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;
398 }
399 }
400 }
401 faceNECount += 2 * local_NE1 * local_NE0;
402 }
403 totalNECount += 2 * NE1 * NE0;
404 }
405 if (local_NE0 > 0)
406 {
407 /* ** elements on boundary 001 (x1=0): - Left */
408
409 if (e_offset0 == 0)
410 {
411 #pragma omp parallel for private(i1,i2,k,node0)
412 for (i2 = 0; i2 < local_NE2; i2++)
413 {
414 for (i1 = 0; i1 < local_NE1; i1++)
415 {
416
417 index_t res, n0, n2, n4, n6;
418 k = 2 * (i1 + local_NE1 * i2) + faceNECount;
419 node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2);
420 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
421 out->FaceElements->Id[k] = res;
422 out->FaceElements->Tag[k] = LEFTTAG;
423 out->FaceElements->Owner[k] = myRank;
424 out->FaceElements->Id[k + 1] = res + 1;
425 out->FaceElements->Tag[k + 1] = LEFTTAG;
426 out->FaceElements->Owner[k + 1] = myRank;
427
428 n0 = node0;
429 n2 = node0 + Nstride1;
430 n4 = node0 + Nstride2;
431 n6 = node0 + Nstride1 + Nstride2;
432
433 if ((global_adjustment + 0 + i1 + i2) % 2 == 0)
434 {
435 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
436 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;
437 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;
438
439 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
440 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;
441 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;
442 }
443 else
444 {
445 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
446
447 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
448 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;
449 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2;
450
451 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;
452 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;
453 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;
454 }
455 }
456 }
457 faceNECount += 2 * local_NE1 * local_NE2;
458 }
459 totalNECount += 2 * NE1 * NE2;
460 /* ** elements on boundary 002 (x1=1): - Right */
461 if (local_NE0 + e_offset0 == NE0)
462 {
463 #pragma omp parallel for private(i1,i2,k,node0)
464 for (i2 = 0; i2 < local_NE2; i2++)
465 {
466 for (i1 = 0; i1 < local_NE1; i1++)
467 {
468 index_t res, n1, n3, n5, n7;
469 k = 2 * (i1 + local_NE1 * i2) + faceNECount;
470
471 node0 =
472 Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) +
473 Nstride2 * N_PER_E * (i2 + e_offset2);
474 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
475 out->FaceElements->Id[k] = res;
476 out->FaceElements->Tag[k] = RIGHTTAG;
477 out->FaceElements->Owner[k] = myRank;
478 out->FaceElements->Id[k + 1] = res + 1;
479 out->FaceElements->Tag[k + 1] = RIGHTTAG;
480 out->FaceElements->Owner[k + 1] = myRank;
481
482 n1 = node0 + Nstride0;
483 n3 = node0 + Nstride0 + Nstride1;
484 n5 = node0 + Nstride0 + Nstride2;
485 n7 = node0 + Nstride0 + Nstride1 + Nstride2;
486 if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0)
487 {
488 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
489 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;
490 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
491
492 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n3;
493 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
494 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5;
495 }
496 else
497 {
498 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
499
500 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
501 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7;
502 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
503
504 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
505 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n3;
506 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n7;
507 }
508 }
509 }
510 faceNECount += 2 * local_NE1 * local_NE2;
511 }
512 totalNECount += 2 * NE1 * NE2;
513 }
514 if (local_NE1 > 0)
515 {
516 /* ** elements on boundary 010 (x2=0): -Front */
517 if (e_offset1 == 0)
518 {
519 #pragma omp parallel for private(i0,i2,k,node0)
520 for (i2 = 0; i2 < local_NE2; i2++)
521 {
522 for (i0 = 0; i0 < local_NE0; i0++)
523 {
524 index_t res, n0, n1, n4, n5;
525 k = 2 * (i0 + local_NE0 * i2) + faceNECount;
526 node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2);
527 res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount;
528 out->FaceElements->Id[k] = res;
529 out->FaceElements->Tag[k] = FRONTTAG;
530 out->FaceElements->Owner[k] = myRank;
531 out->FaceElements->Id[k + 1] = res + 1;
532 out->FaceElements->Tag[k + 1] = FRONTTAG;
533 out->FaceElements->Owner[k + 1] = myRank;
534
535 n0 = node0;
536 n1 = node0 + Nstride0;
537 n4 = node0 + Nstride2;
538 n5 = node0 + Nstride0 + Nstride2;
539
540 if ((global_adjustment + i0 + 0 + i2) % 2 == 0)
541 {
542 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
543 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;
544 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
545
546 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
547 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;
548 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;
549
550 }
551 else
552 {
553 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
554
555 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
556 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;
557 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4;
558
559 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
560 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;
561 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;
562
563 }
564 }
565 }
566 faceNECount += 2 * local_NE0 * local_NE2;
567 }
568 totalNECount += 2 * NE0 * NE2;
569 /* ** elements on boundary 020 (x2=1): - Back */
570 if (local_NE1 + e_offset1 == NE1)
571 {
572 #pragma omp parallel for private(i0,i2,k,node0)
573 for (i2 = 0; i2 < local_NE2; i2++)
574 {
575 for (i0 = 0; i0 < local_NE0; i0++)
576 {
577 index_t res, n2, n6, n7, n3;
578 k = 2 * (i0 + local_NE0 * i2) + faceNECount;
579 node0 =
580 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1) +
581 Nstride2 * N_PER_E * (i2 + e_offset2);
582 res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount;
583 out->FaceElements->Id[k] = res;
584 out->FaceElements->Tag[k] = BACKTAG;
585 out->FaceElements->Owner[k] = myRank;
586 out->FaceElements->Id[k + 1] = res + 1;
587 out->FaceElements->Tag[k + 1] = BACKTAG;
588 out->FaceElements->Owner[k + 1] = myRank;
589
590 n2 = node0 + Nstride1;
591 n6 = node0 + Nstride1 + Nstride2;
592 n7 = node0 + Nstride0 + Nstride1 + Nstride2;
593 n3 = node0 + Nstride0 + Nstride1;
594
595 if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0)
596 {
597 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
598 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;
599 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n3;
600
601 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n6;
602 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
603 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
604
605 }
606 else
607 {
608 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
609 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
610 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;
611 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;
612
613 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2;
614 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
615 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
616
617 }
618 }
619 }
620 faceNECount += 2 * local_NE0 * local_NE2;
621 }
622 totalNECount += 2 * NE0 * NE2;
623 }
624 }
625 if (Dudley_noError())
626 {
627 /* add tag names */
628 Dudley_Mesh_addTagMap(out, "top", TOPTAG);
629 Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG);
630 Dudley_Mesh_addTagMap(out, "left", LEFTTAG);
631 Dudley_Mesh_addTagMap(out, "right", RIGHTTAG);
632 Dudley_Mesh_addTagMap(out, "front", FRONTTAG);
633 Dudley_Mesh_addTagMap(out, "back", BACKTAG);
634 }
635 /* prepare mesh for further calculatuions: */
636 if (Dudley_noError())
637 {
638 Dudley_Mesh_resolveNodeIds(out);
639 }
640 if (Dudley_noError())
641 {
642 Dudley_Mesh_prepare(out, optimize);
643 }
644
645 if (!Dudley_noError())
646 {
647 Dudley_Mesh_free(out);
648 }
649 /* free up memory */
650 Esys_MPIInfo_free(mpi_info);
651
652 return out;
653 }

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