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

Annotation of /trunk/finley/src/Mesh_hex8.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_hex8.c
File MIME type: text/plain
File size: 13882 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 first order elements (Hex8) 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     /* Copyrights by ACcESS Australia 2003/04 */
12     /* Author: gross@access.edu.au */
13     /* Version: $Id$ */
14    
15     /**************************************************************/
16    
17     #include "Common.h"
18     #include "Finley.h"
19     #include "Mesh.h"
20     #include "RectangularMesh.h"
21    
22     /**************************************************************/
23    
24     Finley_Mesh* Finley_RectangularMesh_Hex8(int* numElements,double* Length,int* periodic, int order,int useElementsOnFace) {
25     int N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,node0,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;
26     Finley_Mesh* out;
27     char name[50];
28     double time0=Finley_timer();
29     NE0=MAX(1,numElements[0]);
30     NE1=MAX(1,numElements[1]);
31     NE2=MAX(1,numElements[2]);
32     N0=NE0+1;
33     N1=NE1+1;
34     N2=NE2+1;
35    
36     NFaceElements=0;
37     if (!periodic[0]) {
38     NDOF0=N0;
39     NFaceElements+=2*NE1*NE2;
40     } else {
41     NDOF0=N0-1;
42     }
43     if (!periodic[1]) {
44     NDOF1=N1;
45     NFaceElements+=2*NE0*NE2;
46     } else {
47     NDOF1=N1-1;
48     }
49     if (!periodic[2]) {
50     NDOF2=N2;
51     NFaceElements+=2*NE0*NE1;
52     } else {
53     NDOF2=N2-1;
54     }
55    
56     /* allocate mesh: */
57    
58     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
59     out=Finley_Mesh_alloc(name,3,order);
60     if (Finley_ErrorCode!=NO_ERROR) return NULL;
61    
62     out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
63     if (useElementsOnFace) {
64     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
65     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
66     } else {
67     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
68     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
69     }
70     out->Points=Finley_ElementFile_alloc(Point1,out->order);
71     if (Finley_ErrorCode!=NO_ERROR) {
72     Finley_Mesh_dealloc(out);
73     return NULL;
74     }
75    
76    
77     /* allocate tables: */
78    
79     Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
80     Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
81     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
82     if (Finley_ErrorCode!=NO_ERROR) {
83     Finley_Mesh_dealloc(out);
84     return NULL;
85     }
86    
87     /* set nodes: */
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     /* tags for the faces: */
104     if (!periodic[2]) {
105     for (i1=0;i1<N1;i1++) {
106     for (i0=0;i0<N0;i0++) {
107     out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;
108     out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;
109     }
110     }
111     }
112     if (!periodic[1]) {
113     for (i2=0;i2<N2;i2++) {
114     for (i0=0;i0<N0;i0++) {
115     out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;
116     out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;
117     }
118     }
119     }
120     if (!periodic[0]) {
121     for (i2=0;i2<N2;i2++) {
122     for (i1=0;i1<N1;i1++) {
123     out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;
124     out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;
125     }
126     }
127     }
128    
129     /* set the elements: */
130    
131     #pragma omp parallel for private(i0,i1,i2,k,node0)
132     for (i2=0;i2<NE2;i2++) {
133     for (i1=0;i1<NE1;i1++) {
134     for (i0=0;i0<NE0;i0++) {
135     k=i0+NE0*i1+NE0*NE1*i2;
136     node0=i0+i1*N0+N0*N1*i2;
137    
138     out->Elements->Id[k]=k;
139     out->Elements->Tag[k]=0;
140     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
141    
142     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
143     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
144     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
145     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
146     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
147     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
148     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
149     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
150    
151     }
152     }
153     }
154     out->Elements->numColors=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0)+1;
155    
156     /* face elements: */
157    
158     if (useElementsOnFace) {
159     NUMNODES=8;
160     } else {
161     NUMNODES=4;
162     }
163     totalNECount=NE0*NE1*NE2;
164     faceNECount=0;
165    
166     /* these are the quadrilateral elements on boundary 1 (x3=0): */
167     if (!periodic[2]) {
168     /* ** elements on boundary 100 (x3=0): */
169    
170     #pragma omp parallel for private(i0,i1,k,node0)
171     for (i1=0;i1<NE1;i1++) {
172     for (i0=0;i0<NE0;i0++) {
173     k=i0+NE0*i1+faceNECount;
174     node0=i0+i1*N0;
175    
176     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
177     out->FaceElements->Tag[k]=100;
178     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
179    
180     if (useElementsOnFace) {
181     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
182     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
183     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
184     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
185     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
186     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
187     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
188     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
189     } else {
190     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
191     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
192     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
193     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
194     }
195     }
196     }
197     totalNECount+=NE1*NE0;
198     faceNECount+=NE1*NE0;
199    
200     /* ** elements on boundary 200 (x3=1) */
201    
202     #pragma omp parallel for private(i0,i1,k,node0)
203     for (i1=0;i1<NE1;i1++) {
204     for (i0=0;i0<NE0;i0++) {
205     k=i0+NE0*i1+faceNECount;
206     node0=i0+i1*N0+N0*N1*(NE2-1);
207    
208     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
209     out->FaceElements->Tag[k]=200;
210     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
211    
212     if (useElementsOnFace) {
213     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
214     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
215     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
216     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
217     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
218     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
219     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
220     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
221     } else {
222     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
223     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
224     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
225     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
226     }
227    
228    
229     }
230     }
231     totalNECount+=NE1*NE0;
232     faceNECount+=NE1*NE0;
233     }
234     if (!periodic[0]) {
235     /* ** elements on boundary 001 (x1=0): */
236    
237     #pragma omp parallel for private(i1,i2,k,node0)
238     for (i2=0;i2<NE2;i2++) {
239     for (i1=0;i1<NE1;i1++) {
240     k=i1+NE1*i2+faceNECount;
241     node0=i1*N0+N0*N1*i2;
242    
243     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
244     out->FaceElements->Tag[k]=1;
245     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
246    
247     if (useElementsOnFace) {
248     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
249     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
250     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
251     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
252     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
253     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
254     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
255     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
256     } else {
257     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
258     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
259     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
260     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
261     }
262     }
263     }
264     totalNECount+=NE1*NE2;
265     faceNECount+=NE1*NE2;
266    
267     /* ** elements on boundary 002 (x1=1): */
268    
269     #pragma omp parallel for private(i1,i2,k,node0)
270     for (i2=0;i2<NE2;i2++) {
271     for (i1=0;i1<NE1;i1++) {
272     k=i1+NE1*i2+faceNECount;
273     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
274    
275     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
276     out->FaceElements->Tag[k]=2;
277     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
278    
279     if (useElementsOnFace) {
280     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
281     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
282     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
283     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
284     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
285     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
286     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
287     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
288     } else {
289     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
290     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
291     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
292     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
293     }
294     }
295     }
296     totalNECount+=NE1*NE2;
297     faceNECount+=NE1*NE2;
298     }
299     if (!periodic[1]) {
300     /* ** elements on boundary 010 (x2=0): */
301    
302     #pragma omp parallel for private(i0,i2,k,node0)
303     for (i2=0;i2<NE2;i2++) {
304     for (i0=0;i0<NE0;i0++) {
305     k=i0+NE0*i2+faceNECount;
306     node0=i0+N0*N1*i2;
307    
308     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
309     out->FaceElements->Tag[k]=10;
310     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
311    
312     if (useElementsOnFace) {
313     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
314     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
315     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
316     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
317     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
318     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
319     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
320     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
321     } else {
322     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
323     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
324     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
325     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
326     }
327     }
328     }
329     totalNECount+=NE0*NE2;
330     faceNECount+=NE0*NE2;
331    
332     /* ** elements on boundary 020 (x2=1): */
333    
334     #pragma omp parallel for private(i0,i2,k,node0)
335     for (i2=0;i2<NE2;i2++) {
336     for (i0=0;i0<NE0;i0++) {
337     k=i0+NE0*i2+faceNECount;
338     node0=i0+(NE1-1)*N0+N0*N1*i2;
339    
340     out->FaceElements->Tag[k]=20;
341     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
342     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
343    
344     if (useElementsOnFace) {
345     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
346     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
347     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
348     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
349     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
350     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
351     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
352     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
353     } else {
354     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
355     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
356     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
357     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
358     }
359    
360     }
361     }
362     totalNECount+=NE0*NE2;
363     faceNECount+=NE0*NE2;
364     }
365     out->FaceElements->numColors=24;
366    
367     /* face elements done: */
368    
369     /* condense the nodes: */
370    
371     Finley_Mesh_resolveNodeIds(out);
372    
373     /* prepare mesh for further calculatuions:*/
374     Finley_Mesh_prepare(out) ;
375    
376     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
377    
378     if (Finley_ErrorCode!=NO_ERROR) {
379     Finley_Mesh_dealloc(out);
380     return NULL;
381     }
382     return out;
383     }
384    
385     /*
386     * $Log$
387     * Revision 1.1 2004/10/26 06:53:57 jgs
388     * Initial revision
389     *
390     * Revision 1.1.1.1 2004/06/24 04:00:40 johng
391     * Initial version of eys using boost-python.
392     *
393     *
394     */
395    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26