/[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 6079 - (show annotations)
Mon Mar 21 12:22:38 2016 UTC (2 years, 10 months ago) by caltinay
File size: 25450 byte(s)
Big commit - making dudley much more like finley to make it more
managable. Fixed quite a few issues that had been fixed in finley.
Disposed of all ReducedNode/ReducedDOF entities that dudley never supported.
Compiles and passes tests.

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

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