/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_tet4.cpp
ViewVC logotype

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

Parent Directory Parent Directory | Revision Log Revision Log


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 * http://www.uq.edu.au
6 *
7 * Primary Business: Queensland, Australia
8 * Licensed under the Open Software License version 3.0
9 * http://www.opensource.org/licenses/osl-3.0.php
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