/[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 3117 - (show annotations)
Mon Aug 30 02:10:26 2010 UTC (8 years, 9 months ago) by jfenwick
File MIME type: text/plain
File size: 25228 byte(s)
Face elts fixed

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 Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,
29 double* Length,
30 index_t order,
31 index_t reduced_order,
32 bool_t optimize)
33 {
34 #define N_PER_E 1
35 #define DIM 3
36 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;
37 dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN;
38 index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2;
39 Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL;
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 refElements= Dudley_ReferenceElementSet_alloc(Tet4,order,reduced_order);
81 refFaceElements=Dudley_ReferenceElementSet_alloc(Tri3, order, reduced_order);
82 refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);
83
84
85 if ( Dudley_noError()) {
86
87 Dudley_Mesh_setPoints(out,Dudley_ElementFile_alloc(refPoints, mpi_info));
88 Dudley_Mesh_setFaceElements(out,Dudley_ElementFile_alloc(refFaceElements, mpi_info));
89 Dudley_Mesh_setElements(out,Dudley_ElementFile_alloc(refElements, mpi_info));
90
91 /* work out the largest dimension */
92 if (N2==MAX3(N0,N1,N2)) {
93 Nstride0=1;
94 Nstride1=N0;
95 Nstride2=N0*N1;
96 local_NE0=NE0;
97 e_offset0=0;
98 local_NE1=NE1;
99 e_offset1=0;
100 Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
101 } else if (N1==MAX3(N0,N1,N2)) {
102 Nstride0=N2;
103 Nstride1=N0*N2;
104 Nstride2=1;
105 local_NE0=NE0;
106 e_offset0=0;
107 Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
108 local_NE2=NE2;
109 e_offset2=0;
110 } else {
111 Nstride0=N1*N2;
112 Nstride1=1;
113 Nstride2=N1;
114 Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
115 local_NE1=NE1;
116 e_offset1=0;
117 local_NE2=NE2;
118 e_offset2=0;
119 }
120 offset0=e_offset0*N_PER_E;
121 offset1=e_offset1*N_PER_E;
122 offset2=e_offset2*N_PER_E;
123 local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
124 local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
125 local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
126
127 /* get the number of surface elements */
128
129 NFaceElements=0;
130 if (local_NE2>0) {
131 NDOF2=N2;
132 if (offset2==0) NFaceElements+=2*local_NE1*local_NE0; // each face is split
133 if (local_NE2+e_offset2 == NE2) NFaceElements+=2*local_NE1*local_NE0;
134 } else {
135 NDOF2=N2-1;
136 }
137
138 if (local_NE0>0) {
139 NDOF0=N0;
140 if (e_offset0 == 0) NFaceElements+=2*local_NE1*local_NE2;
141 if (local_NE0+e_offset0 == NE0) NFaceElements+=2*local_NE1*local_NE2;
142 } else {
143 NDOF0=N0-1;
144 }
145
146 if (local_NE1>0) {
147 NDOF1=N1;
148 if (e_offset1 == 0) NFaceElements+=2*local_NE0*local_NE2;
149 if (local_NE1+e_offset1 == NE1) NFaceElements+=2*local_NE0*local_NE2;
150 } else {
151 NDOF1=N1-1;
152 }
153 }
154
155
156 /* allocate tables: */
157 if (Dudley_noError()) {
158
159 Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
160 /* we split the rectangular prism this code used to produce into 5 tetrahedrons */
161 Dudley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2*5);
162 /* each border face will be split in half */
163 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
164 }
165
166 if (Dudley_noError()) {
167 /* create nodes */
168
169 #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
170 for (i2=0;i2<local_N2;i2++) {
171 for (i1=0;i1<local_N1;i1++) {
172 for (i0=0;i0<local_N0;i0++) {
173 k=i0+local_N0*i1+local_N0*local_N1*i2;
174 global_i0=i0+offset0;
175 global_i1=i1+offset1;
176 global_i2=i2+offset2;
177 out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
178 out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
179 out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
180 out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
181 out->Nodes->Tag[k]=0;
182 out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
183 +Nstride1*(global_i1%NDOF1)
184 +Nstride2*(global_i2%NDOF2);
185 /* printf("N=%d: %f,%f,%f\n", k, out->Nodes->Coordinates[INDEX2(0,k,DIM)],
186 out->Nodes->Coordinates[INDEX2(1,k,DIM)],
187 out->Nodes->Coordinates[INDEX2(2,k,DIM)]);*/
188 }
189 }
190 }
191 // printf("Now listing elements\n");
192 /* set the elements: */
193
194 int global_adjustment=0; // If we are not the only rank we may need to shift our pattern to match neighbours
195
196 NN=out->Elements->numNodes;
197 #pragma omp parallel for private(i0,i1,i2,k,node0)
198 for (i2=0;i2<local_NE2;i2++) {
199 for (i1=0;i1<local_NE1;i1++) {
200 for (i0=0;i0<local_NE0;i0++) {
201
202 k=5*(i0+local_NE0*i1+local_NE0*local_NE1*i2);
203 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
204
205 index_t res=5*((i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2));
206 for (int j=0;j<5;++j)
207 {
208 out->Elements->Id[k+j]=res+j;
209 out->Elements->Tag[k+j]=0;
210 out->Elements->Owner[k+j]=myRank;
211 }
212
213 index_t v[8];
214 /* in non-rotated orientation the points are numbered as follows:
215 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
216 if ((global_adjustment+i0+i1+i2)%2==0)
217 {
218 // printf("Type0 %d, %d, %d\n",i0,i1,i2);
219 v[0]=node0;
220 v[1]=node0+Nstride0;
221 v[2]=node0+Nstride1;
222 v[3]=node0+Nstride1+Nstride0;
223 v[4]=node0+Nstride2;
224 v[5]=node0+Nstride0+Nstride2;
225 v[6]=node0+Nstride1+Nstride2;
226 v[7]=node0+Nstride2+Nstride1+Nstride0;
227 }
228 else
229 {
230 // printf("Type1 %d, %d, %d\n",i0,i1,i2);
231 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
232
233 v[0]=node0+Nstride1; // node 0 ends up in position 2
234 v[2]=node0+Nstride1+Nstride2; // node 2 ends up in position 6
235 v[6]=node0+Nstride2; // node 6 ends up in position 4
236 v[4]=node0; // node 4 ends up in position 0
237
238 v[1]=node0+Nstride0+Nstride1; // node 1 -> pos 3
239 v[3]=node0+Nstride2+Nstride1+Nstride0; // node 3-> pos 7
240 v[7]=node0+Nstride0+Nstride2; // node 7 -> pos 5
241 v[5]=node0+Nstride0; // node 5 -> pos 1
242 }
243
244 // for (int z=0;z<8;++z)
245 // {
246 // printf("z[%d]=%d\n", z, v[z]);
247 //
248 // }
249
250
251 // index_t a=node0, b=node0+Nstride0, c=node0+Nstride1+Nstride0, d=node0+Nstride1;
252 // index_t e=node0+Nstride2, f=node0+Nstride2+Nstride0, g=node0+Nstride2+Nstride1+Nstride0,
253 // h=node0+Nstride2+Nstride1;
254 //
255 //
256 //
257 // a=0, b=1, c=3, d=2
258 // e=4, f=5, g=7, h=6
259
260 // elements nodes are numbered: centre, x, y, z
261
262 out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];
263 out->Elements->Nodes[INDEX2(1,k,NN)] =v[5];
264 out->Elements->Nodes[INDEX2(2,k,NN)] =v[6];
265 out->Elements->Nodes[INDEX2(3,k,NN)] =v[0];
266
267 out->Elements->Nodes[INDEX2(0,k+1,NN)] =v[7];
268 out->Elements->Nodes[INDEX2(1,k+1,NN)] =v[6];
269 out->Elements->Nodes[INDEX2(2,k+1,NN)] =v[5];
270 out->Elements->Nodes[INDEX2(3,k+1,NN)] =v[3];
271
272 out->Elements->Nodes[INDEX2(0,k+2,NN)] =v[2];
273 out->Elements->Nodes[INDEX2(1,k+2,NN)] =v[3];
274 out->Elements->Nodes[INDEX2(2,k+2,NN)] =v[0];
275 out->Elements->Nodes[INDEX2(3,k+2,NN)] =v[6];
276
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 // for (int z=0;z<5;++z)
291 // {
292 // printf("E %d:",z);
293 // for (int q=0;q<4;++q)
294 // {
295 // index_t id=out->Elements->Nodes[INDEX2(q,z,NN)];
296 // printf(" %d = %f, %f, %f\n", id, out->Nodes->Coordinates[INDEX2(0,id,DIM)],
297 // out->Nodes->Coordinates[INDEX2(1,id,DIM)], out->Nodes->Coordinates[INDEX2(2,id,DIM)]);
298 // }
299 // printf("\n");}
300
301
302 /*
303 for (int j=0;j<4;++j)
304 {
305
306 printf("Elt %d",j);
307 for (int m=0;m<4;++m)
308 {
309 printf(" %d",out->Elements->Nodes[INDEX2(m,k+j,NN)]);
310 }
311 printf("\n");
312 }*/
313
314
315
316 }
317
318 }
319 } /* end for */
320 /* face elements */
321 // printf("Starting face elements\n");
322 NN=out->FaceElements->numNodes;
323 totalNECount=5*NE0*NE1*NE2;
324 faceNECount=0;
325 //printf("Top\n");
326 /* these are the quadrilateral elements on boundary 1 (x3=0): */
327 if (local_NE2>0)
328 {
329 /* ** elements on boundary 100 (x3=0): */
330 if (e_offset2==0)
331 {
332 #pragma omp parallel for private(i0,i1,k,node0)
333 for (i1=0;i1<local_NE1;i1++)
334 {
335 for (i0=0;i0<local_NE0;i0++)
336 {
337 k=2*(i0+local_NE0*i1)+faceNECount;
338 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
339
340 index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
341 out->FaceElements->Id[k]=res;
342 out->FaceElements->Tag[k]=TOPTAG;
343 out->FaceElements->Owner[k]=myRank;
344 out->FaceElements->Id[k+1]=res+1;
345 out->FaceElements->Tag[k+1]=TOPTAG;
346 out->FaceElements->Owner[k+1]=myRank;
347 // in the element generation this face is a,b,c,d which is split by a-c
348
349 index_t n4=node0+Nstride2;
350 index_t n5=node0+Nstride0+Nstride2;
351 index_t n6=node0+Nstride1+Nstride2;
352 index_t n7=node0+Nstride0+Nstride1+Nstride2;
353
354 if ((global_adjustment+i0+i1)%2==0)
355 {
356 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
357 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
358 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;
359
360 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n5;
361 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
362 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
363
364
365 }
366 else
367 {
368 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n4;
369 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n5;
370 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;
371
372 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;
373 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
374 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n6;
375
376
377
378 }
379 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
380 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);
381 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
382 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
383
384 }
385 }
386 // printf("\n");
387 faceNECount+=2*local_NE1*local_NE0;
388 }
389 totalNECount+=2*NE1*NE0;
390
391 //printf("Bottom\n");
392 /* ** elements on boundary 200 (x3=1) */
393 if (local_NE2+e_offset2 == NE2) {
394 #pragma omp parallel for private(i0,i1,k,node0)
395 for (i1=0;i1<local_NE1;i1++) {
396 for (i0=0;i0<local_NE0;i0++) {
397
398 k=2*(i0+local_NE0*i1)+faceNECount;
399 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
400
401 index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
402 out->FaceElements->Id[k]=res;
403 out->FaceElements->Tag[k]=BOTTOMTAG;
404 out->FaceElements->Owner[k]=myRank;
405 out->FaceElements->Id[k+1]=res+1;
406 out->FaceElements->Tag[k+1]=BOTTOMTAG;
407 out->FaceElements->Owner[k+1]=myRank;
408
409 index_t n0=node0;
410 index_t n1=node0+Nstride0;
411 index_t n2=node0+Nstride1;
412 index_t n3=node0+Nstride1+Nstride0;
413
414 if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
415 {
416 // out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
417 // out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
418 // out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
419 //
420 // out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
421 // out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[1];
422 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
423
424 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
425 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;
426 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;
427
428 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
429 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n1;
430 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
431
432
433 }
434 else
435 {
436
437 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
438 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;
439 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;
440
441 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
442 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;
443 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
444 }
445
446 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
447 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);
448 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
449 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
450
451 }
452 }
453 /*printf("\n");*/
454 faceNECount+=2*local_NE1*local_NE0;
455 }
456 totalNECount+=2*NE1*NE0;
457 }
458 //printf("Left\n");
459 if (local_NE0>0) {
460 /* ** elements on boundary 001 (x1=0): */
461
462 if (e_offset0 == 0) {
463 #pragma omp parallel for private(i1,i2,k,node0)
464 for (i2=0;i2<local_NE2;i2++) {
465 for (i1=0;i1<local_NE1;i1++) {
466
467 k=2*(i1+local_NE1*i2)+faceNECount;
468
469 //printf("%d, %d ",k,k+1);
470 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
471 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
472 out->FaceElements->Id[k]=res;
473 out->FaceElements->Tag[k]=LEFTTAG;
474 out->FaceElements->Owner[k]=myRank;
475 out->FaceElements->Id[k+1]=res+1;
476 out->FaceElements->Tag[k+1]=LEFTTAG;
477 out->FaceElements->Owner[k+1]=myRank;
478
479 index_t n0=node0;
480 index_t n2=node0+Nstride1;
481 index_t n4=node0+Nstride2;
482 index_t n6=node0+Nstride1+Nstride2;
483
484 if ((global_adjustment+0+i1+i2)%2==0)
485 {
486 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
487 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4;
488 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n6;
489
490 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
491 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;
492 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
493 }
494 else
495 {
496 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
497
498 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
499 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n4;
500 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n2;
501
502 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n4;
503 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n6;
504 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n2;
505 }
506 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
507 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);
508 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
509 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
510 }
511 }
512 faceNECount+=2*local_NE1*local_NE2;
513 }
514 /*printf("\n");*/
515 //printf("Right\n");
516 totalNECount+=2*NE1*NE2;
517 /* ** elements on boundary 002 (x1=1): */
518 if (local_NE0+e_offset0 == NE0) {
519 #pragma omp parallel for private(i1,i2,k,node0)
520 for (i2=0;i2<local_NE2;i2++) {
521 for (i1=0;i1<local_NE1;i1++) {
522 k=2*(i1+local_NE1*i2)+faceNECount;
523
524 //printf("%d, %d ",k,k+1);
525 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
526 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
527 out->FaceElements->Id[k]=res;
528 out->FaceElements->Tag[k]=RIGHTTAG;
529 out->FaceElements->Owner[k]=myRank;
530 out->FaceElements->Id[k+1]=res+1;
531 out->FaceElements->Tag[k+1]=RIGHTTAG;
532 out->FaceElements->Owner[k+1]=myRank;
533
534 index_t n1=node0+Nstride0;
535 index_t n3=node0+Nstride0+Nstride1;
536 index_t n5=node0+Nstride0+Nstride2;
537 index_t n7=node0+Nstride0+Nstride1+Nstride2;
538 if ((global_adjustment+local_NE0-1+i1+i2)%2==0)
539 {
540 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1;
541 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n3;
542 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
543
544 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n3;
545 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
546 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n5;
547 }
548 else
549 {
550 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
551
552 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n1;
553 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n7;
554 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
555
556 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
557 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n3;
558 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n7;
559 }
560 // printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
561 // out->FaceElements->Nodes[INDEX2(2,k,NN)]);
562 // printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
563 // out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
564
565 }
566 }
567 faceNECount+=2*local_NE1*local_NE2;
568 }
569 totalNECount+=2*NE1*NE2;
570 }
571 //printf("\n");
572 //printf("Front\n");
573 if (local_NE1>0) {
574 /* ** elements on boundary 010 (x2=0): */
575 if (e_offset1 == 0) {
576 #pragma omp parallel for private(i0,i2,k,node0)
577 for (i2=0;i2<local_NE2;i2++) {
578 for (i0=0;i0<local_NE0;i0++) {
579 k=2*(i0+local_NE0*i2)+faceNECount;
580 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
581 //printf("%d, %d ",k,k+1);
582 index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;
583 out->FaceElements->Id[k]=res;
584 out->FaceElements->Tag[k]=FRONTTAG;
585 out->FaceElements->Owner[k]=myRank;
586 out->FaceElements->Id[k+1]=res+1;
587 out->FaceElements->Tag[k+1]=FRONTTAG;
588 out->FaceElements->Owner[k+1]=myRank;
589
590 index_t n0=node0;
591 index_t n1=node0+Nstride0;
592 index_t n4=node0+Nstride2;
593 index_t n5=node0+Nstride0+Nstride2;
594
595 if ((global_adjustment+i0+0+i2)%2==0)
596 {
597 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
598 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;
599 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n5;
600
601 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n0;
602 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5;
603 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;
604
605
606 }
607 else
608 {
609 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
610
611 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n0;
612 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n1;
613 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n4;
614
615 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n1;
616 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n5;
617 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n4;
618
619 }
620 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
621 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
622 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
623 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/
624 }
625 }
626 faceNECount+=2*local_NE0*local_NE2;
627 }
628 // printf("\n");
629 totalNECount+=2*NE0*NE2;
630 //printf("Back\n");
631 /* ** elements on boundary 020 (x2=1): */
632 if (local_NE1+e_offset1 == NE1) {
633 #pragma omp parallel for private(i0,i2,k,node0)
634 for (i2=0;i2<local_NE2;i2++) {
635 for (i0=0;i0<local_NE0;i0++) {
636 k=2*(i0+local_NE0*i2)+faceNECount;
637 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
638 // printf("%d, %d ",k,k+1);
639 index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;
640 out->FaceElements->Id[k]=res;
641 out->FaceElements->Tag[k]=BACKTAG;
642 out->FaceElements->Owner[k]=myRank;
643 out->FaceElements->Id[k+1]=res+1;
644 out->FaceElements->Tag[k+1]=BACKTAG;
645 out->FaceElements->Owner[k+1]=myRank;
646
647 index_t n2=node0+Nstride1;
648 index_t n6=node0+Nstride1+Nstride2;
649 index_t n7=node0+Nstride0+Nstride1+Nstride2;
650 index_t n3=node0+Nstride0+Nstride1;
651
652 if ((global_adjustment+i0+local_NE1-1+i2)%2==0)
653 {
654 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2;
655 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6;
656 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n3;
657
658 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n6;
659 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
660 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
661
662
663 }
664 else
665 {
666 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
667 out->FaceElements->Nodes[INDEX2(0,k,NN)]=n2;
668 out->FaceElements->Nodes[INDEX2(1,k,NN)]=n6;
669 out->FaceElements->Nodes[INDEX2(2,k,NN)]=n7;
670
671 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=n2;
672 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=n7;
673 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=n3;
674
675 }
676 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
677 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
678 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
679 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/
680 }
681 }
682 faceNECount+=2*local_NE0*local_NE2;
683 }
684 totalNECount+=2*NE0*NE2;
685 }
686 }
687 // printf("\n");
688 if (Dudley_noError()) {
689 /* add tag names */
690 Dudley_Mesh_addTagMap(out,"top", TOPTAG);
691 Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
692 Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
693 Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
694 Dudley_Mesh_addTagMap(out,"front", FRONTTAG);
695 Dudley_Mesh_addTagMap(out,"back", BACKTAG);
696 }
697 /* prepare mesh for further calculatuions:*/
698 if (Dudley_noError()) {
699 Dudley_Mesh_resolveNodeIds(out);
700 }
701 if (Dudley_noError()) {
702 Dudley_Mesh_prepare(out, optimize);
703 }
704
705 if (!Dudley_noError()) {
706 Dudley_Mesh_free(out);
707 }
708 /* free up memory */
709 Dudley_ReferenceElementSet_dealloc(refPoints);
710 Dudley_ReferenceElementSet_dealloc(refFaceElements);
711 Dudley_ReferenceElementSet_dealloc(refElements);
712 Paso_MPIInfo_free( mpi_info );
713
714 return out;
715 }
716
717

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26