/[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 3220 - (show annotations)
Wed Sep 29 00:33:16 2010 UTC (8 years, 7 months ago) by jfenwick
File MIME type: text/plain
File size: 24897 byte(s)
Removing references to Quadrature.? and ShapeFns

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26