/[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 3114 - (show annotations)
Fri Aug 27 05:26:25 2010 UTC (8 years, 10 months ago) by jfenwick
File MIME type: text/plain
File size: 25128 byte(s)
It doesn't pass all tests but this is major progress

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 top face (clockwise= 0,1,3,2), bottom face (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+Nstride2;
234 v[2]=node0;
235 v[6]=node0+Nstride1;
236 v[4]=node0+Nstride1+Nstride2;
237
238 v[1]=node0+Nstride0+Nstride2;
239 v[3]=node0+Nstride0;
240 v[7]=node0+Nstride1+Nstride0;
241 v[5]=node0+Nstride2+Nstride1+Nstride0;
242 }
243 // index_t a=node0, b=node0+Nstride0, c=node0+Nstride1+Nstride0, d=node0+Nstride1;
244 // index_t e=node0+Nstride2, f=node0+Nstride2+Nstride0, g=node0+Nstride2+Nstride1+Nstride0,
245 // h=node0+Nstride2+Nstride1;
246 //
247 //
248 //
249 // a=0, b=1, c=3, d=2
250 // e=4, f=5, g=7, h=6
251
252 // elements nodes are numbered: centre, x, y, z
253
254 out->Elements->Nodes[INDEX2(0,k,NN)] =v[4];
255 out->Elements->Nodes[INDEX2(1,k,NN)] =v[5];
256 out->Elements->Nodes[INDEX2(2,k,NN)] =v[6];
257 out->Elements->Nodes[INDEX2(3,k,NN)] =v[0];
258
259 out->Elements->Nodes[INDEX2(0,k+1,NN)] =v[7];
260 out->Elements->Nodes[INDEX2(1,k+1,NN)] =v[6];
261 out->Elements->Nodes[INDEX2(2,k+1,NN)] =v[5];
262 out->Elements->Nodes[INDEX2(3,k+1,NN)] =v[3];
263
264 out->Elements->Nodes[INDEX2(0,k+2,NN)] =v[2];
265 out->Elements->Nodes[INDEX2(1,k+2,NN)] =v[3];
266 out->Elements->Nodes[INDEX2(2,k+2,NN)] =v[0];
267 out->Elements->Nodes[INDEX2(3,k+2,NN)] =v[6];
268
269 out->Elements->Nodes[INDEX2(0,k+3,NN)] =v[1];
270 out->Elements->Nodes[INDEX2(1,k+3,NN)] =v[0];
271 out->Elements->Nodes[INDEX2(2,k+3,NN)] =v[3];
272 out->Elements->Nodes[INDEX2(3,k+3,NN)] =v[5];
273
274 // I can't work out where the center is for this one
275 out->Elements->Nodes[INDEX2(0,k+4,NN)] =v[5];
276 out->Elements->Nodes[INDEX2(1,k+4,NN)] =v[0];
277 out->Elements->Nodes[INDEX2(2,k+4,NN)] =v[6];
278 out->Elements->Nodes[INDEX2(3,k+4,NN)] =v[3];
279
280
281
282 for (int j=0;j<4;++j)
283 {
284
285 printf("Elt %d",j);
286 for (int m=0;m<4;++m)
287 {
288 printf(" %d",out->Elements->Nodes[INDEX2(m,k+j,NN)]);
289 }
290 printf("\n");
291 }
292
293
294
295 }
296
297 }
298 } /* end for */
299 /* face elements */
300 // printf("Starting face elements\n");
301 NN=out->FaceElements->numNodes;
302 totalNECount=5*NE0*NE1*NE2;
303 faceNECount=0;
304 printf("Top\n");
305 /* these are the quadrilateral elements on boundary 1 (x3=0): */
306 if (local_NE2>0)
307 {
308 /* ** elements on boundary 100 (x3=0): */
309 if (e_offset2==0)
310 {
311 #pragma omp parallel for private(i0,i1,k,node0)
312 for (i1=0;i1<local_NE1;i1++)
313 {
314 for (i0=0;i0<local_NE0;i0++)
315 {
316 k=2*(i0+local_NE0*i1)+faceNECount;
317 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
318
319 index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
320 out->FaceElements->Id[k]=res;
321 out->FaceElements->Tag[k]=TOPTAG;
322 out->FaceElements->Owner[k]=myRank;
323 out->FaceElements->Id[k+1]=res+1;
324 out->FaceElements->Tag[k+1]=TOPTAG;
325 out->FaceElements->Owner[k+1]=myRank;
326 // in the element generation this face is a,b,c,d which is split by a-c
327
328 index_t v[4];
329 if ((global_adjustment+i0+i1)%2==0)
330 {
331 v[0]=node0+Nstride1+Nstride2; // node 6
332 v[1]=node0+Nstride2+Nstride1+Nstride0; // node 7
333 v[2]=node0+Nstride2; // node 4
334 v[3]=node0+Nstride0+Nstride2; // node 5
335
336 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
337 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[1];
338 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[3];
339
340 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
341 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
342 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[2];
343
344 }
345 else
346 {
347 v[0]=node0+Nstride2; // node4
348 v[1]=node0+Nstride0+Nstride2; // node5
349 v[2]=node0; // node0
350 v[3]=node0+Nstride0; // node1
351
352 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
353 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[1];
354 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
355
356 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
357 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
358 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
359 }
360 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
361 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
362 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
363 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
364
365 }
366 }
367 // printf("\n");
368 faceNECount+=2*local_NE1*local_NE0;
369 }
370 totalNECount+=2*NE1*NE0;
371
372 printf("Bottom\n");
373 /* ** elements on boundary 200 (x3=1) */
374 if (local_NE2+e_offset2 == NE2) {
375 #pragma omp parallel for private(i0,i1,k,node0)
376 for (i1=0;i1<local_NE1;i1++) {
377 for (i0=0;i0<local_NE0;i0++) {
378
379 k=2*(i0+local_NE0*i1)+faceNECount;
380 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
381
382 index_t res=2*((i0+e_offset0)+NE0*(i1+e_offset1))+totalNECount;
383 out->FaceElements->Id[k]=res;
384 out->FaceElements->Tag[k]=BOTTOMTAG;
385 out->FaceElements->Owner[k]=myRank;
386 out->FaceElements->Id[k+1]=res+1;
387 out->FaceElements->Tag[k+1]=BOTTOMTAG;
388 out->FaceElements->Owner[k+1]=myRank;
389
390 index_t v[4];
391 if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
392 {
393 v[0]=node0; // node 0
394 v[1]=node0+Nstride0; // node 1
395 v[2]=node0+Nstride1; // node 2
396 v[3]=node0+Nstride1+Nstride0; // node 3
397
398 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
399 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
400 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
401
402 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
403 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[1];
404 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
405
406
407
408 }
409 else
410 {
411 v[0]=node0+Nstride2; // 4
412 v[1]=node0+Nstride0+Nstride2; //5
413 v[2]=node0; // 0
414 v[3]=node0+Nstride0; // 1
415
416
417 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
418 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
419 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
420
421 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
422 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
423 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
424 }
425
426 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
427 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
428 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
429 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
430
431 }
432 }
433 /*printf("\n");*/
434 faceNECount+=2*local_NE1*local_NE0;
435 }
436 totalNECount+=2*NE1*NE0;
437 }
438 printf("Left\n");
439 if (local_NE0>0) {
440 /* ** elements on boundary 001 (x1=0): */
441
442 if (e_offset0 == 0) {
443 #pragma omp parallel for private(i1,i2,k,node0)
444 for (i2=0;i2<local_NE2;i2++) {
445 for (i1=0;i1<local_NE1;i1++) {
446
447 k=2*(i1+local_NE1*i2)+faceNECount;
448
449 //printf("%d, %d ",k,k+1);
450 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
451 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
452 out->FaceElements->Id[k]=res;
453 out->FaceElements->Tag[k]=LEFTTAG;
454 out->FaceElements->Owner[k]=myRank;
455 out->FaceElements->Id[k+1]=res+1;
456 out->FaceElements->Tag[k+1]=LEFTTAG;
457 out->FaceElements->Owner[k+1]=myRank;
458
459 index_t v[4];
460 if ((global_adjustment+0+i1+i2)%2==0)
461 {
462 v[0]=node0+Nstride1+Nstride2; // 6
463 v[1]=node0+Nstride2; // 4
464 v[2]=node0+Nstride1; // 2
465 v[3]=node0;
466 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
467 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
468 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
469
470 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
471 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
472 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
473 }
474 else
475 {
476 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
477 v[0]=node0+Nstride1; // 2
478 v[1]=node0+Nstride1+Nstride2; // 6
479 v[2]=node0; // 0
480 v[3]=node0+Nstride2; // 4
481
482 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
483 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
484 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
485
486 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
487 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
488 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
489 }
490 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
491 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
492 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
493 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
494 }
495 }
496 faceNECount+=2*local_NE1*local_NE2;
497 }
498 /*printf("\n");*/
499 printf("Right\n");
500 totalNECount+=2*NE1*NE2;
501 /* ** elements on boundary 002 (x1=1): */
502 if (local_NE0+e_offset0 == NE0) {
503 #pragma omp parallel for private(i1,i2,k,node0)
504 for (i2=0;i2<local_NE2;i2++) {
505 for (i1=0;i1<local_NE1;i1++) {
506 k=2*(i1+local_NE1*i2)+faceNECount;
507
508 //printf("%d, %d ",k,k+1);
509 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
510 index_t res=2*((i1+e_offset1)+NE1*(i2+e_offset2))+totalNECount;
511 out->FaceElements->Id[k]=res;
512 out->FaceElements->Tag[k]=RIGHTTAG;
513 out->FaceElements->Owner[k]=myRank;
514 out->FaceElements->Id[k+1]=res+1;
515 out->FaceElements->Tag[k+1]=RIGHTTAG;
516 out->FaceElements->Owner[k+1]=myRank;
517
518 index_t v[4];
519 if ((global_adjustment+local_NE0-1+i1+i2)%2==0)
520 {
521 v[0]=node0+Nstride0+Nstride2; // 5
522 v[1]=node0+Nstride2+Nstride1+Nstride0; // 7
523 v[2]=node0+Nstride0; // 1
524 v[3]=node0+Nstride1+Nstride0; // 3
525
526
527 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
528 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
529 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
530
531 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
532 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
533 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
534 }
535 else
536 {
537 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
538 v[0]=node0+Nstride2; // 4
539 v[1]=node0+Nstride1+Nstride0; // 3
540 v[2]=node0+Nstride0+Nstride2; // 5
541 v[3]=node0+Nstride0; // 1
542
543 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
544 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
545 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
546
547 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[1];
548 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[2];
549 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[3];
550 }
551 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
552 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
553 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
554 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
555
556 }
557 }
558 faceNECount+=2*local_NE1*local_NE2;
559 }
560 totalNECount+=2*NE1*NE2;
561 }
562 //printf("\n");
563 printf("Front\n");
564 if (local_NE1>0) {
565 /* ** elements on boundary 010 (x2=0): */
566 if (e_offset1 == 0) {
567 #pragma omp parallel for private(i0,i2,k,node0)
568 for (i2=0;i2<local_NE2;i2++) {
569 for (i0=0;i0<local_NE0;i0++) {
570 k=2*(i0+local_NE0*i2)+faceNECount;
571 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
572 //printf("%d, %d ",k,k+1);
573 index_t res=2*((i2+e_offset2)+NE2*(e_offset0+i0))+totalNECount;
574 out->FaceElements->Id[k]=res;
575 out->FaceElements->Tag[k]=FRONTTAG;
576 out->FaceElements->Owner[k]=myRank;
577 out->FaceElements->Id[k+1]=res+1;
578 out->FaceElements->Tag[k+1]=FRONTTAG;
579 out->FaceElements->Owner[k+1]=myRank;
580
581 index_t v[4];
582 if ((global_adjustment+i0+0+i2)%2==0)
583 {
584 v[0]=node0+Nstride2; // 4
585 v[1]=node0+Nstride0+Nstride2; // 5
586 v[2]=node0; // 0
587 v[3]=node0+Nstride0; // 1
588
589
590 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
591 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
592 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
593
594 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
595 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
596 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
597
598
599 }
600 else
601 {
602 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
603 v[0]=node0+Nstride1+Nstride2; // 6
604 v[1]=node0+Nstride2+Nstride1+Nstride0; // 7
605 v[2]=node0+Nstride2; // 4
606 v[3]=node0+Nstride0+Nstride2; // 5
607
608 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
609 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[3];
610 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[2];
611
612 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
613 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
614 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
615
616 }
617 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
618 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
619 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
620 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
621 }
622 }
623 faceNECount+=2*local_NE0*local_NE2;
624 }
625 // printf("\n");
626 totalNECount+=2*NE0*NE2;
627 printf("Back\n");
628 /* ** elements on boundary 020 (x2=1): */
629 if (local_NE1+e_offset1 == NE1) {
630 #pragma omp parallel for private(i0,i2,k,node0)
631 for (i2=0;i2<local_NE2;i2++) {
632 for (i0=0;i0<local_NE0;i0++) {
633 k=2*(i0+local_NE0*i2)+faceNECount;
634 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
635 // printf("%d, %d ",k,k+1);
636 index_t res=2*((i2+e_offset2)+NE2*(i0+e_offset0))+totalNECount;
637 out->FaceElements->Id[k]=res;
638 out->FaceElements->Tag[k]=BACKTAG;
639 out->FaceElements->Owner[k]=myRank;
640 out->FaceElements->Id[k+1]=res+1;
641 out->FaceElements->Tag[k+1]=BACKTAG;
642 out->FaceElements->Owner[k+1]=myRank;
643
644 index_t v[4];
645 if ((global_adjustment+i0+local_NE1-1+i2)%2==0)
646 {
647 v[0]=node0+Nstride2+Nstride1+Nstride0; // 7
648 v[1]=node0+Nstride1+Nstride2; // 6
649 v[2]=node0+Nstride1+Nstride0; // 3
650 v[3]=node0+Nstride1; // 2
651
652
653 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
654 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
655 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[1];
656
657 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[2];
658 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
659 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
660
661
662 }
663 else
664 {
665 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
666 v[0]=node0+Nstride1+Nstride0; // 3
667 v[1]=node0+Nstride1; // 2
668 v[2]=node0+Nstride0; // 1
669 v[3]=node0; // 0
670
671 out->FaceElements->Nodes[INDEX2(0,k,NN)]=v[0];
672 out->FaceElements->Nodes[INDEX2(1,k,NN)]=v[2];
673 out->FaceElements->Nodes[INDEX2(2,k,NN)]=v[3];
674
675 out->FaceElements->Nodes[INDEX2(0,k+1,NN)]=v[0];
676 out->FaceElements->Nodes[INDEX2(1,k+1,NN)]=v[3];
677 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]=v[1];
678
679 }
680 printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
681 out->FaceElements->Nodes[INDEX2(2,k,NN)]);
682 printf("%d: %d, %d, %d\n", k+1,out->FaceElements->Nodes[INDEX2(0,k+1,NN)], out->FaceElements->Nodes[INDEX2(1,k+1,NN)],
683 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);
684 }
685 }
686 faceNECount+=2*local_NE0*local_NE2;
687 }
688 totalNECount+=2*NE0*NE2;
689 }
690 }
691 // printf("\n");
692 if (Dudley_noError()) {
693 /* add tag names */
694 Dudley_Mesh_addTagMap(out,"top", TOPTAG);
695 Dudley_Mesh_addTagMap(out,"bottom", BOTTOMTAG);
696 Dudley_Mesh_addTagMap(out,"left", LEFTTAG);
697 Dudley_Mesh_addTagMap(out,"right", RIGHTTAG);
698 Dudley_Mesh_addTagMap(out,"front", FRONTTAG);
699 Dudley_Mesh_addTagMap(out,"back", BACKTAG);
700 }
701 /* prepare mesh for further calculatuions:*/
702 if (Dudley_noError()) {
703 Dudley_Mesh_resolveNodeIds(out);
704 }
705 if (Dudley_noError()) {
706 Dudley_Mesh_prepare(out, optimize);
707 }
708
709 if (!Dudley_noError()) {
710 Dudley_Mesh_free(out);
711 }
712 /* free up memory */
713 Dudley_ReferenceElementSet_dealloc(refPoints);
714 Dudley_ReferenceElementSet_dealloc(refFaceElements);
715 Dudley_ReferenceElementSet_dealloc(refElements);
716 Paso_MPIInfo_free( mpi_info );
717
718 return out;
719 }
720
721

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26