/[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 123 - (hide annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 22259 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26