/[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 3114 - (hide annotations)
Fri Aug 27 05:26:25 2010 UTC (8 years, 9 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 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 3114 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     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 jgs 82 }
296 jfenwick 3114
297 jgs 82 }
298 jfenwick 3114 } /* end for */
299 ksteube 1312 /* face elements */
300 jfenwick 3114 // printf("Starting face elements\n");
301 ksteube 1312 NN=out->FaceElements->numNodes;
302 jfenwick 3114 totalNECount=5*NE0*NE1*NE2;
303 ksteube 1312 faceNECount=0;
304 jfenwick 3114 printf("Top\n");
305 ksteube 1312 /* these are the quadrilateral elements on boundary 1 (x3=0): */
306 jfenwick 3114 if (local_NE2>0)
307     {
308 ksteube 1312 /* ** elements on boundary 100 (x3=0): */
309 jfenwick 3114 if (e_offset2==0)
310     {
311 ksteube 1312 #pragma omp parallel for private(i0,i1,k,node0)
312 jfenwick 3114 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 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
318 jgs 82
319 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
323 jfenwick 3114 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 ksteube 1312 }
366     }
367 jfenwick 3114 // printf("\n");
368     faceNECount+=2*local_NE1*local_NE0;
369 jgs 82 }
370 jfenwick 3114 totalNECount+=2*NE1*NE0;
371    
372     printf("Bottom\n");
373 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i1)+faceNECount;
380 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
381 jfenwick 3114
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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
386 jfenwick 3114 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 ksteube 1312 }
432     }
433 jfenwick 3114 /*printf("\n");*/
434     faceNECount+=2*local_NE1*local_NE0;
435 jgs 82 }
436 jfenwick 3114 totalNECount+=2*NE1*NE0;
437 jgs 82 }
438 jfenwick 3114 printf("Left\n");
439     if (local_NE0>0) {
440 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
448    
449     //printf("%d, %d ",k,k+1);
450 ksteube 1312 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
451 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
455 jfenwick 3114 out->FaceElements->Id[k+1]=res+1;
456     out->FaceElements->Tag[k+1]=LEFTTAG;
457     out->FaceElements->Owner[k+1]=myRank;
458 ksteube 1312
459 jfenwick 3114 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 ksteube 1312 }
495     }
496 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
497 ksteube 1312 }
498 jfenwick 3114 /*printf("\n");*/
499     printf("Right\n");
500     totalNECount+=2*NE1*NE2;
501 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i1+local_NE1*i2)+faceNECount;
507    
508     //printf("%d, %d ",k,k+1);
509 ksteube 1312 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
510 jfenwick 3114 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
514 jfenwick 3114 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 ksteube 1312 }
557     }
558 jfenwick 3114 faceNECount+=2*local_NE1*local_NE2;
559 jgs 82 }
560 jfenwick 3114 totalNECount+=2*NE1*NE2;
561 jgs 82 }
562 jfenwick 3114 //printf("\n");
563     printf("Front\n");
564     if (local_NE1>0) {
565 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
571 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
572 jfenwick 3114 //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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
577 jfenwick 3114 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 ksteube 1312 }
622     }
623 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
624 ksteube 1312 }
625 jfenwick 3114 // printf("\n");
626     totalNECount+=2*NE0*NE2;
627     printf("Back\n");
628 ksteube 1312 /* ** 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 jfenwick 3114 k=2*(i0+local_NE0*i2)+faceNECount;
634 ksteube 1312 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
635 jfenwick 3114 // 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 ksteube 1312 out->FaceElements->Owner[k]=myRank;
640 jfenwick 3114 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 ksteube 1312 }
685     }
686 jfenwick 3114 faceNECount+=2*local_NE0*local_NE2;
687 ksteube 1312 }
688 jfenwick 3114 totalNECount+=2*NE0*NE2;
689 bcumming 751 }
690 gross 2748 }
691 jfenwick 3114 // printf("\n");
692 jfenwick 3086 if (Dudley_noError()) {
693 ksteube 1312 /* add tag names */
694 jfenwick 3114 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 gross 2748 }
701     /* prepare mesh for further calculatuions:*/
702 jfenwick 3086 if (Dudley_noError()) {
703     Dudley_Mesh_resolveNodeIds(out);
704 gross 2748 }
705 jfenwick 3086 if (Dudley_noError()) {
706     Dudley_Mesh_prepare(out, optimize);
707 bcumming 751 }
708    
709 jfenwick 3086 if (!Dudley_noError()) {
710     Dudley_Mesh_free(out);
711 bcumming 751 }
712 gross 2748 /* free up memory */
713 jfenwick 3086 Dudley_ReferenceElementSet_dealloc(refPoints);
714     Dudley_ReferenceElementSet_dealloc(refFaceElements);
715     Dudley_ReferenceElementSet_dealloc(refElements);
716 ksteube 1312 Paso_MPIInfo_free( mpi_info );
717    
718     return out;
719 bcumming 751 }
720 jfenwick 3114
721    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26