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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3221 - (show annotations)
Wed Sep 29 01:00:21 2010 UTC (8 years, 6 months ago) by jfenwick
File MIME type: text/plain
File size: 20886 byte(s)
Comment stripping

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision
svn:mergeinfo

  ViewVC Help
Powered by ViewVC 1.1.26