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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 3117 - (hide 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 jgs 150
2 ksteube 1312 /*******************************************************
3 ksteube 1811 *
4 jfenwick 2881 * Copyright (c) 2003-2010 by University of Queensland
5 ksteube 1811 * 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 ksteube 1312
14 ksteube 1811
15 jgs 82 /**************************************************************/
16    
17 jfenwick 3086 /* Dudley: generates rectangular meshes */
18 jgs 82
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 jfenwick 3114 #include "TriangularMesh.h"
26 jgs 82
27 jfenwick 3114
28     Dudley_Mesh* Dudley_TriangularMesh_Tet4(dim_t* numElements,
29 ksteube 1312 double* Length,
30 jfenwick 3114 index_t order,
31 ksteube 1312 index_t reduced_order,
32     bool_t optimize)
33 bcumming 751 {
34 ksteube 1312 #define N_PER_E 1
35     #define DIM 3
36 gross 2748 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 jfenwick 3114 Dudley_ReferenceElementSet *refPoints=NULL, *refFaceElements=NULL, *refElements=NULL;
40 jfenwick 3086 Dudley_Mesh* out;
41 ksteube 1312 Paso_MPIInfo *mpi_info = NULL;
42     char name[50];
43 jfenwick 3086 #ifdef Dudley_TRACE
44     double time0=Dudley_timer();
45 jfenwick 1981 #endif
46 bcumming 751
47 jfenwick 3114 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 ksteube 1312 /* get MPI information */
58     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
59 jfenwick 3086 if (! Dudley_noError()) {
60 ksteube 1312 return NULL;
61 bcumming 751 }
62 ksteube 1312 myRank=mpi_info->rank;
63 bcumming 751
64 ksteube 1312 /* set up the global dimensions of the mesh */
65 bcumming 751
66 jgs 82 NE0=MAX(1,numElements[0]);
67     NE1=MAX(1,numElements[1]);
68     NE2=MAX(1,numElements[2]);
69 ksteube 1312 N0=N_PER_E*NE0+1;
70     N1=N_PER_E*NE1+1;
71     N2=N_PER_E*NE2+1;
72 jgs 82
73 ksteube 1312 /* allocate mesh: */
74 jfenwick 3114 sprintf(name,"Triangular %d x %d x %d (x 5) mesh",N0,N1,N2);
75 jfenwick 3086 out=Dudley_Mesh_alloc(name,DIM, mpi_info);
76     if (! Dudley_noError()) {
77 ksteube 1312 Paso_MPIInfo_free( mpi_info );
78     return NULL;
79 jgs 153 }
80 jfenwick 3114 refElements= Dudley_ReferenceElementSet_alloc(Tet4,order,reduced_order);
81     refFaceElements=Dudley_ReferenceElementSet_alloc(Tri3, order, reduced_order);
82 jfenwick 3086 refPoints=Dudley_ReferenceElementSet_alloc(Point1, order, reduced_order);
83 gross 2748
84    
85 jfenwick 3086 if ( Dudley_noError()) {
86 gross 2748
87 jfenwick 3086 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 ksteube 1312
91 gross 2748 /* 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 bcumming 751
127 gross 2748 /* get the number of surface elements */
128 jgs 82
129 gross 2748 NFaceElements=0;
130 jfenwick 3114 if (local_NE2>0) {
131 gross 2748 NDOF2=N2;
132 jfenwick 3114 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 gross 2748 } else {
135     NDOF2=N2-1;
136     }
137 ksteube 1312
138 jfenwick 3114 if (local_NE0>0) {
139 gross 2748 NDOF0=N0;
140 jfenwick 3114 if (e_offset0 == 0) NFaceElements+=2*local_NE1*local_NE2;
141     if (local_NE0+e_offset0 == NE0) NFaceElements+=2*local_NE1*local_NE2;
142 gross 2748 } else {
143     NDOF0=N0-1;
144     }
145 jfenwick 3114
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 jgs 82 }
154 gross 2748
155    
156 ksteube 1312 /* allocate tables: */
157 jfenwick 3086 if (Dudley_noError()) {
158 jgs 82
159 jfenwick 3086 Dudley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
160 jfenwick 3114 /* 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 jfenwick 3086 Dudley_ElementFile_allocTable(out->FaceElements,NFaceElements);
164 gross 2748 }
165    
166 jfenwick 3086 if (Dudley_noError()) {
167 ksteube 1312 /* create nodes */
168 jgs 82
169 ksteube 1312 #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 jfenwick 3117 /* 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 jgs 82 }
189     }
190     }
191 jfenwick 3114 // printf("Now listing elements\n");
192 ksteube 1312 /* set the elements: */
193 jfenwick 3114
194     int global_adjustment=0; // If we are not the only rank we may need to shift our pattern to match neighbours
195    
196 ksteube 1312 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 jfenwick 3114 k=5*(i0+local_NE0*i1+local_NE0*local_NE1*i2);
203 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
204 jgs 82
205 jfenwick 3114 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 jfenwick 3117 The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
216 jfenwick 3114 if ((global_adjustment+i0+i1+i2)%2==0)
217     {
218 jfenwick 3115 // printf("Type0 %d, %d, %d\n",i0,i1,i2);
219 jfenwick 3114 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 jfenwick 3115 // printf("Type1 %d, %d, %d\n",i0,i1,i2);
231 jfenwick 3114 // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
232    
233 jfenwick 3117 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 jfenwick 3114
238 jfenwick 3117 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 jfenwick 3114 }
243 jfenwick 3117
244     // for (int z=0;z<8;++z)
245     // {
246     // printf("z[%d]=%d\n", z, v[z]);
247     //
248     // }
249    
250    
251 jfenwick 3114 // 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 jfenwick 3117
278 jfenwick 3114 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 jfenwick 3117 // 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 jfenwick 3115 /*
303 jfenwick 3114 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 jfenwick 3115 }*/
313 jfenwick 3114
314    
315    
316 jgs 82 }
317 jfenwick 3114
318 jgs 82 }
319 jfenwick 3114 } /* end for */
320 ksteube 1312 /* face elements */
321 jfenwick 3114 // printf("Starting face elements\n");
322 ksteube 1312 NN=out->FaceElements->numNodes;
323 jfenwick 3114 totalNECount=5*NE0*NE1*NE2;
324 ksteube 1312 faceNECount=0;
325 jfenwick 3115 //printf("Top\n");
326 ksteube 1312 /* these are the quadrilateral elements on boundary 1 (x3=0): */
327 jfenwick 3114 if (local_NE2>0)
328     {
329 ksteube 1312 /* ** elements on boundary 100 (x3=0): */
330 jfenwick 3114 if (e_offset2==0)
331     {
332 ksteube 1312 #pragma omp parallel for private(i0,i1,k,node0)
333 jfenwick 3114 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 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
339 jgs 82
340 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
344 jfenwick 3114 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 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+i0+i1)%2==0)
355     {
356 jfenwick 3117 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 jfenwick 3114
360 jfenwick 3117 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 jfenwick 3114
364    
365     }
366     else
367     {
368 jfenwick 3117 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 jfenwick 3114
372 jfenwick 3117 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 jfenwick 3114
376 jfenwick 3117
377    
378 jfenwick 3114 }
379 jfenwick 3115 // 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 jfenwick 3114
384 ksteube 1312 }
385     }
386 jfenwick 3114 // printf("\n");
387     faceNECount+=2*local_NE1*local_NE0;
388 jgs 82 }
389 jfenwick 3114 totalNECount+=2*NE1*NE0;
390    
391 jfenwick 3115 //printf("Bottom\n");
392 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i1)+faceNECount;
399 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
400 jfenwick 3114
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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
405 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
406     out->FaceElements->Tag[k+1]=BOTTOMTAG;
407     out->FaceElements->Owner[k+1]=myRank;
408    
409 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+i0+i1+local_NE2-1)%2==0)
415     {
416 jfenwick 3117 // 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 jfenwick 3114
424 jfenwick 3117 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 jfenwick 3114
428 jfenwick 3117 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 jfenwick 3114
432    
433     }
434     else
435     {
436    
437 jfenwick 3117 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 jfenwick 3114
441 jfenwick 3117 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 jfenwick 3114 }
445    
446 jfenwick 3115 // 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 jfenwick 3114
451 ksteube 1312 }
452     }
453 jfenwick 3114 /*printf("\n");*/
454     faceNECount+=2*local_NE1*local_NE0;
455 jgs 82 }
456 jfenwick 3114 totalNECount+=2*NE1*NE0;
457 jgs 82 }
458 jfenwick 3115 //printf("Left\n");
459 jfenwick 3114 if (local_NE0>0) {
460 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
468    
469     //printf("%d, %d ",k,k+1);
470 ksteube 1312 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
471 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
475 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
476     out->FaceElements->Tag[k+1]=LEFTTAG;
477     out->FaceElements->Owner[k+1]=myRank;
478 ksteube 1312
479 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+0+i1+i2)%2==0)
485     {
486 jfenwick 3117 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 jfenwick 3114
490 jfenwick 3117 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 jfenwick 3114 }
494     else
495     {
496     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
497    
498 jfenwick 3117 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 jfenwick 3114
502 jfenwick 3117 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 jfenwick 3114 }
506 jfenwick 3115 // 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 ksteube 1312 }
511     }
512 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
513 ksteube 1312 }
514 jfenwick 3114 /*printf("\n");*/
515 jfenwick 3115 //printf("Right\n");
516 jfenwick 3114 totalNECount+=2*NE1*NE2;
517 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
523    
524     //printf("%d, %d ",k,k+1);
525 ksteube 1312 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
526 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
530 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
531     out->FaceElements->Tag[k+1]=RIGHTTAG;
532     out->FaceElements->Owner[k+1]=myRank;
533    
534 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+local_NE0-1+i1+i2)%2==0)
539     {
540 jfenwick 3117 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 jfenwick 3114
544 jfenwick 3117 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 jfenwick 3114 }
548     else
549     {
550     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
551    
552 jfenwick 3117 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 jfenwick 3114
556 jfenwick 3117 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 jfenwick 3114 }
560 jfenwick 3115 // 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 jfenwick 3114
565 ksteube 1312 }
566     }
567 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
568 jgs 82 }
569 jfenwick 3114 totalNECount+=2*NE1*NE2;
570 jgs 82 }
571 jfenwick 3114 //printf("\n");
572 jfenwick 3115 //printf("Front\n");
573 jfenwick 3114 if (local_NE1>0) {
574 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
580 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
581 jfenwick 3114 //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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
586 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
587     out->FaceElements->Tag[k+1]=FRONTTAG;
588     out->FaceElements->Owner[k+1]=myRank;
589    
590 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+i0+0+i2)%2==0)
596     {
597 jfenwick 3117 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 jfenwick 3114
601 jfenwick 3117 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 jfenwick 3114
605    
606     }
607     else
608     {
609     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
610    
611 jfenwick 3117 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 jfenwick 3114
615 jfenwick 3117 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 jfenwick 3114
619     }
620 jfenwick 3115 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
621 jfenwick 3114 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 jfenwick 3115 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/
624 ksteube 1312 }
625     }
626 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
627 ksteube 1312 }
628 jfenwick 3114 // printf("\n");
629     totalNECount+=2*NE0*NE2;
630 jfenwick 3115 //printf("Back\n");
631 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
637 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
638 jfenwick 3114 // 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
643 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
644     out->FaceElements->Tag[k+1]=BACKTAG;
645     out->FaceElements->Owner[k+1]=myRank;
646    
647 jfenwick 3117 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 jfenwick 3114 if ((global_adjustment+i0+local_NE1-1+i2)%2==0)
653     {
654 jfenwick 3117 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 jfenwick 3114
658 jfenwick 3117 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 jfenwick 3114
662    
663     }
664     else
665     {
666     // this form is rotated around the 0,2,4,6 face clockwise 90 degrees
667 jfenwick 3117 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 jfenwick 3114
671 jfenwick 3117 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 jfenwick 3114
675     }
676 jfenwick 3115 /*printf("%d: %d, %d, %d\n", k,out->FaceElements->Nodes[INDEX2(0,k,NN)], out->FaceElements->Nodes[INDEX2(1,k,NN)],
677 jfenwick 3114 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 jfenwick 3115 out->FaceElements->Nodes[INDEX2(2,k+1,NN)]);*/
680 ksteube 1312 }
681     }
682 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
683 ksteube 1312 }
684 jfenwick 3114 totalNECount+=2*NE0*NE2;
685 bcumming 751 }
686 gross 2748 }
687 jfenwick 3114 // printf("\n");
688 jfenwick 3086 if (Dudley_noError()) {
689 ksteube 1312 /* add tag names */
690 jfenwick 3114 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 gross 2748 }
697     /* prepare mesh for further calculatuions:*/
698 jfenwick 3086 if (Dudley_noError()) {
699     Dudley_Mesh_resolveNodeIds(out);
700 gross 2748 }
701 jfenwick 3086 if (Dudley_noError()) {
702     Dudley_Mesh_prepare(out, optimize);
703 bcumming 751 }
704    
705 jfenwick 3086 if (!Dudley_noError()) {
706     Dudley_Mesh_free(out);
707 bcumming 751 }
708 gross 2748 /* free up memory */
709 jfenwick 3086 Dudley_ReferenceElementSet_dealloc(refPoints);
710     Dudley_ReferenceElementSet_dealloc(refFaceElements);
711     Dudley_ReferenceElementSet_dealloc(refElements);
712 ksteube 1312 Paso_MPIInfo_free( mpi_info );
713    
714     return out;
715 bcumming 751 }
716 jfenwick 3114
717    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26