/[escript]/trunk/finley/src/Mesh_hex20.c
ViewVC logotype

Annotation of /trunk/finley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 21953 byte(s)
Initial revision

1 jgs 82 /**************************************************************/
2    
3     /* Finley: generates rectangular meshes */
4    
5     /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with second order elements (Hex20) in the brick */
6     /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
7     /* integration scheme. */
8    
9    
10     /**************************************************************/
11    
12     /* Copyrights by ACcESS Australia 2003/04 */
13     /* Author: gross@access.edu.au */
14     /* Version: $Id$ */
15    
16     /**************************************************************/
17    
18     #include "Common.h"
19     #include "Finley.h"
20     #include "Mesh.h"
21     #include "RectangularMesh.h"
22    
23     /**************************************************************/
24    
25     Finley_Mesh* Finley_RectangularMesh_Hex20(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
26     int N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,node0,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;
27     Finley_Mesh* out;
28     char name[50];
29     double time0=Finley_timer();
30    
31     NE0=MAX(1,numElements[0]);
32     NE1=MAX(1,numElements[1]);
33     NE2=MAX(1,numElements[2]);
34     N0=2*NE0+1;
35     N1=2*NE1+1;
36     N2=2*NE2+1;
37    
38     NFaceElements=0;
39     if (!periodic[0]) {
40     NDOF0=N0;
41     NFaceElements+=2*NE1*NE2;
42     } else {
43     NDOF0=N0-1;
44     }
45     if (!periodic[1]) {
46     NDOF1=N1;
47     NFaceElements+=2*NE0*NE2;
48     } else {
49     NDOF1=N1-1;
50     }
51     if (!periodic[2]) {
52     NDOF2=N2;
53     NFaceElements+=2*NE0*NE1;
54     } else {
55     NDOF2=N2-1;
56     }
57    
58     /* allocate mesh: */
59    
60     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
61     out=Finley_Mesh_alloc(name,3,order);
62     if (Finley_ErrorCode!=NO_ERROR) return NULL;
63    
64     out->Elements=Finley_ElementFile_alloc(Hex20,out->order);
65     if (useElementsOnFace) {
66     out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order);
67     out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order);
68     } else {
69     out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order);
70     out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order);
71     }
72     out->Points=Finley_ElementFile_alloc(Point1,out->order);
73     if (Finley_ErrorCode!=NO_ERROR) {
74     Finley_Mesh_dealloc(out);
75     return NULL;
76     }
77    
78    
79     /* allocate tables: */
80    
81     Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
82     Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
83     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
84     if (Finley_ErrorCode!=NO_ERROR) {
85     Finley_Mesh_dealloc(out);
86     return NULL;
87     }
88    
89     #pragma omp parallel for private(i0,i1,i2,k)
90     for (i2=0;i2<N2;i2++) {
91     for (i1=0;i1<N1;i1++) {
92     for (i0=0;i0<N0;i0++) {
93     k=i0+N0*i1+N0*N1*i2;
94     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
95     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
96     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
97     out->Nodes->Id[k]=k;
98     out->Nodes->Tag[k]=0;
99     out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1) +N0*N1*(i2%NDOF2);
100     }
101     }
102     }
103    
104     /* tags for the faces: */
105     if (!periodic[2]) {
106     for (i1=0;i1<N1;i1++) {
107     for (i0=0;i0<N0;i0++) {
108     out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;
109     out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;
110     }
111     }
112     }
113     if (!periodic[1]) {
114     for (i2=0;i2<N2;i2++) {
115     for (i0=0;i0<N0;i0++) {
116     out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;
117     out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;
118     }
119     }
120     }
121     if (!periodic[0]) {
122     for (i2=0;i2<N2;i2++) {
123     for (i1=0;i1<N1;i1++) {
124     out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;
125     out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;
126     }
127     }
128     }
129    
130     /* set the elements: */
131    
132     #pragma omp parallel for private(i0,i1,i2,k,node0)
133     for (i2=0;i2<NE2;i2++) {
134     for (i1=0;i1<NE1;i1++) {
135     for (i0=0;i0<NE0;i0++) {
136     k=i0+NE0*i1+NE0*NE1*i2;
137     node0=2*i0+2*i1*N0+2*N0*N1*i2;
138    
139     out->Elements->Id[k]=k;
140     out->Elements->Tag[k]=0;
141     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
142    
143     out->Elements->Nodes[INDEX2(0,k,20)]=node0;
144     out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
145     out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;
146     out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;
147     out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;
148     out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;
149     out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;
150     out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;
151     out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
152     out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;
153     out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;
154     out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;
155     out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;
156     out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;
157     out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;
158     out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;
159     out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;
160     out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;
161     out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;
162     out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;
163     }
164     }
165     }
166     out->Elements->numColors=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0)+1;
167    
168     /* face elements: */
169    
170     if (useElementsOnFace) {
171     NUMNODES=20;
172     } else {
173     NUMNODES=8;
174     }
175     totalNECount=NE0*NE1*NE2;
176     faceNECount=0;
177    
178     /* these are the quadrilateral elements on boundary 1 (x3=0): */
179    
180     if (!periodic[2]) {
181     /* ** elements on boundary 100 (x3=0): */
182     #pragma omp parallel for private(i0,i1,k,node0)
183     for (i1=0;i1<NE1;i1++) {
184     for (i0=0;i0<NE0;i0++) {
185     k=i0+NE0*i1+faceNECount;
186     node0=2*i0+2*i1*N0;
187    
188     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
189     out->FaceElements->Tag[k]=100;
190     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
191    
192     if (useElementsOnFace) {
193     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
194     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
195     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
196     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
197     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;
198     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;
199     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
200     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;
201     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;
202     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;
203     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;
204     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;
205     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
206     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;
207     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
208     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;
209     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;
210     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
211     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;
212     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;
213     } else {
214     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
215     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
216     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
217     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
218     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
219     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
220     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
221     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
222     }
223     }
224     }
225     totalNECount+=NE1*NE0;
226     faceNECount+=NE1*NE0;
227    
228     /* ** elements on boundary 200 (x3=1) */
229     #pragma omp parallel for private(i0,i1,k,node0)
230     for (i1=0;i1<NE1;i1++) {
231     for (i0=0;i0<NE0;i0++) {
232     k=i0+NE0*i1+faceNECount;
233     node0=2*i0+2*i1*N0+2*N0*N1*(NE2-1);
234    
235     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
236     out->FaceElements->Tag[k]=200;
237     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
238    
239     if (useElementsOnFace) {
240     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
241     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
242     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
243     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
244    
245     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
246     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;
247     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;
248     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;
249    
250     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;
251     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;
252     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
253     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;
254    
255     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
256     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;
257     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
258     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;
259    
260     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;
261     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;
262     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;
263     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;
264    
265     } else {
266     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
267     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
268     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
269     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
270     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;
271     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;
272     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
273     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;
274     }
275     }
276     }
277     totalNECount+=NE1*NE0;
278     faceNECount+=NE1*NE0;
279     }
280     if (!periodic[0]) {
281     /* ** elements on boundary 001 (x1=0): */
282    
283     #pragma omp parallel for private(i1,i2,k,node0)
284     for (i2=0;i2<NE2;i2++) {
285     for (i1=0;i1<NE1;i1++) {
286     k=i1+NE1*i2+faceNECount;
287     node0=2*i1*N0+2*N0*N1*i2;
288    
289     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
290     out->FaceElements->Tag[k]=1;
291     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
292    
293     if (useElementsOnFace) {
294     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
295     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
296     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
297     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
298    
299     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;
300     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;
301     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
302     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;
303    
304     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;
305     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;
306     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;
307     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;
308    
309     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
310     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;
311     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
312     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;
313    
314     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;
315     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;
316     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;
317     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;
318    
319     } else {
320     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
321     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
322     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
323     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
324     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
325     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;
326     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;
327     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
328     }
329     }
330     }
331     totalNECount+=NE1*NE2;
332     faceNECount+=NE1*NE2;
333    
334     /* ** elements on boundary 002 (x1=1): */
335    
336     #pragma omp parallel for private(i1,i2,k,node0)
337     for (i2=0;i2<NE2;i2++) {
338     for (i1=0;i1<NE1;i1++) {
339     k=i1+NE1*i2+faceNECount;
340     node0=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;
341    
342     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
343     out->FaceElements->Tag[k]=2;
344     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
345    
346     if (useElementsOnFace) {
347     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
348     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
349     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
350     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
351    
352     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
353     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;
354     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;
355     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;
356    
357     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;
358     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;
359     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;
360     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;
361    
362     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
363     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;
364     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
365     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;
366    
367     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;
368     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;
369     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;
370     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;
371    
372     } else {
373     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
374     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
375     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
376     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
377     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
378     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;
379     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;
380     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;
381     }
382    
383     }
384     }
385     totalNECount+=NE1*NE2;
386     faceNECount+=NE1*NE2;
387     }
388     if (!periodic[1]) {
389     /* ** elements on boundary 010 (x2=0): */
390    
391     #pragma omp parallel for private(i0,i2,k,node0)
392     for (i2=0;i2<NE2;i2++) {
393     for (i0=0;i0<NE0;i0++) {
394     k=i0+NE0*i2+faceNECount;
395     node0=2*i0+2*N0*N1*i2;
396    
397     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
398     out->FaceElements->Tag[k]=10;
399     out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
400    
401     if (useElementsOnFace) {
402     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
403     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
404     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
405     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
406    
407     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;
408     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;
409     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;
410     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;
411    
412     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;
413     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;
414     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;
415     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;
416    
417     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
418     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;
419     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;
420     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;
421    
422     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;
423     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;
424     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
425     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;
426    
427     } else {
428     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
429     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
430     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
431     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
432     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
433     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;
434     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;
435     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;
436     }
437     }
438     }
439     totalNECount+=NE0*NE2;
440     faceNECount+=NE0*NE2;
441    
442     /* ** elements on boundary 020 (x2=1): */
443    
444     #pragma omp parallel for private(i0,i2,k,node0)
445     for (i2=0;i2<NE2;i2++) {
446     for (i0=0;i0<NE0;i0++) {
447     k=i0+NE0*i2+faceNECount;
448     node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;
449    
450     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
451     out->FaceElements->Tag[k]=20;
452     out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
453    
454     if (useElementsOnFace) {
455     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
456     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
457     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
458     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
459    
460     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
461     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;
462     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;
463     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;
464    
465     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;
466     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
467     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;
468     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;
469    
470     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
471     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;
472     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;
473     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;
474    
475     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;
476     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;
477     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;
478     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;
479     } else {
480     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
481     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
482     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
483     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
484     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;
485     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
486     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;
487     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
488     }
489     }
490     }
491     totalNECount+=NE0*NE2;
492     faceNECount+=NE0*NE2;
493     }
494     out->FaceElements->numColors=24;
495    
496     /* face elements done: */
497    
498     /* condense the nodes: */
499    
500     Finley_Mesh_resolveNodeIds(out);
501    
502     /* prepare mesh for further calculatuions:*/
503     Finley_Mesh_prepare(out) ;
504    
505     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
506    
507     if (Finley_ErrorCode!=NO_ERROR) {
508     Finley_Mesh_dealloc(out);
509     return NULL;
510     }
511     return out;
512     }
513    
514     /*
515     * $Log$
516     * Revision 1.1 2004/10/26 06:53:57 jgs
517     * Initial revision
518     *
519     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
520     * Initial version of eys using boost-python.
521     *
522     *
523     */
524    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26