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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3792 - (hide annotations)
Wed Feb 1 06:16:25 2012 UTC (7 years, 2 months ago) by caltinay
File MIME type: text/plain
File size: 21169 byte(s)
Merged ripley rectangular domain into trunk.

1 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 jgs 82 /**************************************************************/
15    
16 jfenwick 3086 /* Dudley: generates rectangular meshes */
17 jgs 82
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 jfenwick 3114 #include "TriangularMesh.h"
25 jgs 82
26 jfenwick 3118 /* Be careful reading this function. The X? and NStride? are 1,2,3 but the loop vars are 0,1,2 */
27 jfenwick 3224 Dudley_Mesh *Dudley_TriangularMesh_Tet4(dim_t * numElements,
28     double *Length, index_t order, index_t reduced_order, bool_t optimize)
29 bcumming 751 {
30 jfenwick 3224 #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 jfenwick 3227 Esys_MPIInfo *mpi_info = NULL;
39 jfenwick 3224 char name[50];
40     #ifdef Dudley_TRACE
41     double time0 = Dudley_timer();
42     #endif
43 bcumming 751
44 jfenwick 3224 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 jfenwick 3114
51 jfenwick 3224 /* get MPI information */
52 jfenwick 3227 mpi_info = Esys_MPIInfo_alloc(MPI_COMM_WORLD);
53 jfenwick 3224 if (!Dudley_noError())
54     {
55     return NULL;
56     }
57     myRank = mpi_info->rank;
58 jfenwick 3114
59 jfenwick 3224 /* set up the global dimensions of the mesh */
60 jfenwick 3114
61 jfenwick 3224 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 jfenwick 3114
68 jfenwick 3224 /* 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 jfenwick 3227 Esys_MPIInfo_free(mpi_info);
74 jfenwick 3224 return NULL;
75     }
76     if (Dudley_noError())
77     {
78 bcumming 751
79 caltinay 3247 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 bcumming 751
83 jfenwick 3224 /* 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 jfenwick 3227 Esys_MPIInfo_Split(mpi_info, NE2, &local_NE2, &e_offset2);
94 jfenwick 3224 }
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 jfenwick 3227 Esys_MPIInfo_Split(mpi_info, NE1, &local_NE1, &e_offset1);
103 jfenwick 3224 local_NE2 = NE2;
104     e_offset2 = 0;
105     }
106     else
107     {
108     Nstride0 = N1 * N2;
109     Nstride1 = 1;
110     Nstride2 = N1;
111 jfenwick 3227 Esys_MPIInfo_Split(mpi_info, NE0, &local_NE0, &e_offset0);
112 jfenwick 3224 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 jgs 82
124 jfenwick 3224 /* get the number of surface elements */
125 ksteube 1312
126 jfenwick 3224 NFaceElements = 0;
127     if (local_NE2 > 0)
128     {
129     NDOF2 = N2;
130     if (offset2 == 0)
131 caltinay 3254 NFaceElements += 2 * local_NE1 * local_NE0; /* each face is split */
132 jfenwick 3224 if (local_NE2 + e_offset2 == NE2)
133     NFaceElements += 2 * local_NE1 * local_NE0;
134     }
135     else
136     {
137     NDOF2 = N2 - 1;
138     }
139 bcumming 751
140 jfenwick 3224 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 jgs 82
153 jfenwick 3224 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 jgs 82
167 jfenwick 3224 /* allocate tables: */
168     if (Dudley_noError())
169     {
170    
171     Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1 * local_N2);
172 jfenwick 3114 /* we split the rectangular prism this code used to produce into 5 tetrahedrons */
173 jfenwick 3224 Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * local_NE2 * 5);
174 jfenwick 3114 /* each border face will be split in half */
175 jfenwick 3224 Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);
176     }
177 jfenwick 3114
178 jfenwick 3224 if (Dudley_noError())
179     {
180 caltinay 3254 int global_adjustment;
181    
182 jfenwick 3224 /* create nodes */
183 jfenwick 3114
184 jfenwick 3224 #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 jgs 82
207 caltinay 3254 global_adjustment = (offset0 + offset1 + offset2) % 2; /* If we are not the only rank we may need to shift our pattern to match neighbours */
208 jfenwick 3114
209 jfenwick 3224 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 caltinay 3254 index_t res;
218     index_t v[8];
219     int j;
220 jfenwick 3224 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 caltinay 3254 res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2));
226     for (j = 0; j < 5; ++j)
227 jfenwick 3224 {
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 jfenwick 3114 /* in non-rotated orientation the points are numbered as follows:
234 jfenwick 3117 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
235 jfenwick 3224 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 caltinay 3254 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
249 jfenwick 3114
250 caltinay 3254 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 jfenwick 3114
255 caltinay 3254 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 jfenwick 3224 }
260 jfenwick 3117
261 caltinay 3254 /* elements nodes are numbered: centre, x, y, z */
262 jfenwick 3114
263 jfenwick 3224 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 jfenwick 3114
268 jfenwick 3224 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 jfenwick 3114
273 jfenwick 3224 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 jfenwick 3114
278 jfenwick 3224 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 jfenwick 3117
283 caltinay 3254 /* I can't work out where the center is for this one */
284 jfenwick 3224 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 jfenwick 3114
290 jfenwick 3224 }
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 jfenwick 3114 {
302 jfenwick 3224 #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 caltinay 3254 index_t res, n0, n1, n2, n3;
308 jfenwick 3224 k = 2 * (i0 + local_NE0 * i1) + faceNECount;
309     node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);
310 caltinay 3254 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
311 jfenwick 3224 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 jfenwick 3114
318 caltinay 3254 n0 = node0;
319     n1 = node0 + Nstride0;
320     n2 = node0 + Nstride1;
321     n3 = node0 + Nstride0 + Nstride1;
322 jfenwick 3117
323 jfenwick 3224 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 jfenwick 3114
329 jfenwick 3224 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 jfenwick 3114
333 jfenwick 3224 }
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 jfenwick 3114
340 jfenwick 3224 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 jfenwick 3114
344 jfenwick 3224 }
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 jfenwick 3117
359 caltinay 3254 index_t res, n4, n5, n6, n7;
360 jfenwick 3224 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 jfenwick 3117
365 caltinay 3254 res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
366 jfenwick 3224 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 jfenwick 3114
373 caltinay 3254 n4 = node0 + Nstride2;
374     n5 = node0 + Nstride0 + Nstride2;
375     n6 = node0 + Nstride1 + Nstride2;
376     n7 = node0 + Nstride1 + Nstride0 + Nstride2;
377 jfenwick 3114
378 jfenwick 3224 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 jfenwick 3117
384 jfenwick 3224 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 jfenwick 3114
391 jfenwick 3224 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 jfenwick 3114
395 jfenwick 3224 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 jfenwick 3114
409 jfenwick 3224 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 jfenwick 3114
417 caltinay 3254 index_t res, n0, n2, n4, n6;
418 jfenwick 3224 k = 2 * (i1 + local_NE1 * i2) + faceNECount;
419     node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2);
420 caltinay 3254 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
421 jfenwick 3224 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 jfenwick 3114
428 caltinay 3254 n0 = node0;
429     n2 = node0 + Nstride1;
430     n4 = node0 + Nstride2;
431     n6 = node0 + Nstride1 + Nstride2;
432 jfenwick 3114
433 jfenwick 3224 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 jfenwick 3114
439 jfenwick 3224 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 caltinay 3254 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
446 jfenwick 3114
447 jfenwick 3224 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 jfenwick 3114
451 jfenwick 3224 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 caltinay 3254 index_t res, n1, n3, n5, n7;
469 jfenwick 3224 k = 2 * (i1 + local_NE1 * i2) + faceNECount;
470 jfenwick 3114
471 jfenwick 3224 node0 =
472     Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) +
473     Nstride2 * N_PER_E * (i2 + e_offset2);
474 caltinay 3254 res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
475 jfenwick 3224 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 jfenwick 3114
482 caltinay 3254 n1 = node0 + Nstride0;
483     n3 = node0 + Nstride0 + Nstride1;
484     n5 = node0 + Nstride0 + Nstride2;
485     n7 = node0 + Nstride0 + Nstride1 + Nstride2;
486 jfenwick 3224 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 jfenwick 3114
492 jfenwick 3224 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 caltinay 3254 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
499 jfenwick 3117
500 jfenwick 3224 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 jfenwick 3114
504 jfenwick 3224 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 caltinay 3254 index_t res, n0, n1, n4, n5;
525 jfenwick 3224 k = 2 * (i0 + local_NE0 * i2) + faceNECount;
526     node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2);
527 caltinay 3254 res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount;
528 jfenwick 3224 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 jfenwick 3114
535 caltinay 3254 n0 = node0;
536     n1 = node0 + Nstride0;
537     n4 = node0 + Nstride2;
538     n5 = node0 + Nstride0 + Nstride2;
539 jfenwick 3114
540 jfenwick 3224 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 jfenwick 3114
546 jfenwick 3224 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 jfenwick 3114
550 jfenwick 3224 }
551     else
552     {
553 caltinay 3254 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
554 jfenwick 3114
555 jfenwick 3224 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 jfenwick 3114
559 jfenwick 3224 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 jfenwick 3117
563 jfenwick 3224 }
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 caltinay 3254 index_t res, n2, n6, n7, n3;
578 jfenwick 3224 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 caltinay 3254 res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount;
583 jfenwick 3224 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 jfenwick 3114
590 caltinay 3254 n2 = node0 + Nstride1;
591     n6 = node0 + Nstride1 + Nstride2;
592     n7 = node0 + Nstride0 + Nstride1 + Nstride2;
593     n3 = node0 + Nstride0 + Nstride1;
594 jfenwick 3114
595 jfenwick 3224 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 jfenwick 3114
601 jfenwick 3224 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 jfenwick 3114
605 jfenwick 3224 }
606     else
607     {
608 caltinay 3254 /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
609 jfenwick 3224 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 jfenwick 3114
613 jfenwick 3224 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 bcumming 751
617 jfenwick 3224 }
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 gross 2748 /* free up memory */
650 jfenwick 3227 Esys_MPIInfo_free(mpi_info);
651 ksteube 1312
652 jfenwick 3224 return out;
653 bcumming 751 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26