/[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 616 - (hide annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 10 months ago) by elspeth
File MIME type: text/plain
File size: 14782 byte(s)
Copyright added to more source files.

1 jgs 150 /*
2 elspeth 616 ************************************************************
3     * Copyright 2006 by ACcESS MNRF *
4     * *
5     * http://www.access.edu.au *
6     * Primary Business: Queensland, Australia *
7     * Licensed under the Open Software License version 3.0 *
8     * http://www.opensource.org/licenses/osl-3.0.php *
9     * *
10     ************************************************************
11 jgs 150 */
12    
13 jgs 82 /**************************************************************/
14    
15     /* Finley: generates rectangular meshes */
16    
17     /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
18     /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
19     /* integration scheme. */
20    
21     /**************************************************************/
22    
23 jgs 150 /* Author: gross@access.edu.au */
24     /* Version: $Id$ */
25 jgs 82
26     /**************************************************************/
27    
28     #include "RectangularMesh.h"
29    
30     /**************************************************************/
31    
32 jgs 123 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
33 jgs 153 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
34 jgs 123 index_t node0;
35 jgs 82 Finley_Mesh* out;
36     char name[50];
37     double time0=Finley_timer();
38     NE0=MAX(1,numElements[0]);
39     NE1=MAX(1,numElements[1]);
40     NE2=MAX(1,numElements[2]);
41     N0=NE0+1;
42     N1=NE1+1;
43     N2=NE2+1;
44    
45 jgs 153 if (N0<=MIN(N1,N2)) {
46     if (N1 <= N2) {
47     M0=1;
48     M1=N0;
49     M2=N0*N1;
50     } else {
51     M0=1;
52     M2=N0;
53     M1=N0*N2;
54     }
55     } else if (N1<=MIN(N2,N0)) {
56     if (N2 <= N0) {
57     M1=1;
58     M2=N1;
59     M0=N2*N1;
60     } else {
61     M1=1;
62     M0=N1;
63     M2=N1*N0;
64     }
65     } else {
66     if (N0 <= N1) {
67     M2=1;
68     M0=N2;
69     M1=N2*N0;
70     } else {
71     M2=1;
72     M1=N2;
73     M0=N1*N2;
74     }
75     }
76    
77    
78 jgs 82 NFaceElements=0;
79     if (!periodic[0]) {
80     NDOF0=N0;
81     NFaceElements+=2*NE1*NE2;
82     } else {
83     NDOF0=N0-1;
84     }
85     if (!periodic[1]) {
86     NDOF1=N1;
87     NFaceElements+=2*NE0*NE2;
88     } else {
89     NDOF1=N1-1;
90     }
91     if (!periodic[2]) {
92     NDOF2=N2;
93     NFaceElements+=2*NE0*NE1;
94     } else {
95     NDOF2=N2-1;
96     }
97    
98     /* allocate mesh: */
99    
100     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
101     out=Finley_Mesh_alloc(name,3,order);
102 jgs 150 if (! Finley_noError()) return NULL;
103 jgs 82
104     out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
105     if (useElementsOnFace) {
106     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
107     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
108     } else {
109     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
110     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
111     }
112     out->Points=Finley_ElementFile_alloc(Point1,out->order);
113 jgs 150 if (! Finley_noError()) {
114 jgs 82 Finley_Mesh_dealloc(out);
115     return NULL;
116     }
117    
118    
119     /* allocate tables: */
120    
121     Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
122     Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
123     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
124 jgs 150 if (! Finley_noError()) {
125 jgs 82 Finley_Mesh_dealloc(out);
126     return NULL;
127     }
128    
129     /* set nodes: */
130 jgs 153
131     #pragma omp parallel for private(i0,i1,i2,k)
132 jgs 82 for (i2=0;i2<N2;i2++) {
133     for (i1=0;i1<N1;i1++) {
134     for (i0=0;i0<N0;i0++) {
135 jgs 153 k=M0*i0+M1*i1+M2*i2;
136 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
137     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
138     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
139 jgs 153 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
140 jgs 82 out->Nodes->Tag[k]=0;
141 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
142 jgs 82 }
143     }
144     }
145     /* tags for the faces: */
146 jgs 153 /* tags for the faces: */
147 jgs 82 if (!periodic[2]) {
148 jgs 153 for (i1=0;i1<N1;i1++) {
149     for (i0=0;i0<N0;i0++) {
150     out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
151     out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
152     }
153     }
154 jgs 82 }
155     if (!periodic[1]) {
156     for (i2=0;i2<N2;i2++) {
157     for (i0=0;i0<N0;i0++) {
158 jgs 153 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
159     out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
160 jgs 82 }
161     }
162     }
163     if (!periodic[0]) {
164     for (i2=0;i2<N2;i2++) {
165     for (i1=0;i1<N1;i1++) {
166 jgs 153 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
167     out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
168 jgs 82 }
169     }
170     }
171 jgs 153
172 jgs 82 /* set the elements: */
173    
174     #pragma omp parallel for private(i0,i1,i2,k,node0)
175     for (i2=0;i2<NE2;i2++) {
176     for (i1=0;i1<NE1;i1++) {
177     for (i0=0;i0<NE0;i0++) {
178     k=i0+NE0*i1+NE0*NE1*i2;
179     node0=i0+i1*N0+N0*N1*i2;
180    
181     out->Elements->Id[k]=k;
182     out->Elements->Tag[k]=0;
183     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
184    
185     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
186     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
187     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
188     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
189     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
190     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
191     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
192     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
193    
194     }
195     }
196     }
197 jgs 123 out->Elements->minColor=0;
198     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
199 jgs 82
200     /* face elements: */
201    
202     if (useElementsOnFace) {
203     NUMNODES=8;
204     } else {
205     NUMNODES=4;
206     }
207     totalNECount=NE0*NE1*NE2;
208     faceNECount=0;
209    
210     /* these are the quadrilateral elements on boundary 1 (x3=0): */
211     if (!periodic[2]) {
212     /* ** elements on boundary 100 (x3=0): */
213    
214     #pragma omp parallel for private(i0,i1,k,node0)
215     for (i1=0;i1<NE1;i1++) {
216     for (i0=0;i0<NE0;i0++) {
217     k=i0+NE0*i1+faceNECount;
218     node0=i0+i1*N0;
219    
220     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
221     out->FaceElements->Tag[k]=100;
222     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
223    
224     if (useElementsOnFace) {
225     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
226     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
227     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
228     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
229     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
230     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
231     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
232     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
233     } else {
234     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
235     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
236     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
237     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
238     }
239     }
240     }
241     totalNECount+=NE1*NE0;
242     faceNECount+=NE1*NE0;
243    
244     /* ** elements on boundary 200 (x3=1) */
245    
246     #pragma omp parallel for private(i0,i1,k,node0)
247     for (i1=0;i1<NE1;i1++) {
248     for (i0=0;i0<NE0;i0++) {
249     k=i0+NE0*i1+faceNECount;
250     node0=i0+i1*N0+N0*N1*(NE2-1);
251    
252     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
253     out->FaceElements->Tag[k]=200;
254     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
255    
256     if (useElementsOnFace) {
257     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
258     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
259     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
260     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
261     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
262     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
263     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
264     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
265     } else {
266     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
267     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
268     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
269     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
270     }
271    
272    
273     }
274     }
275     totalNECount+=NE1*NE0;
276     faceNECount+=NE1*NE0;
277     }
278     if (!periodic[0]) {
279     /* ** elements on boundary 001 (x1=0): */
280    
281     #pragma omp parallel for private(i1,i2,k,node0)
282     for (i2=0;i2<NE2;i2++) {
283     for (i1=0;i1<NE1;i1++) {
284     k=i1+NE1*i2+faceNECount;
285     node0=i1*N0+N0*N1*i2;
286    
287     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
288     out->FaceElements->Tag[k]=1;
289     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
290    
291     if (useElementsOnFace) {
292     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
293     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
294     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
295     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
296     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
297     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
298     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
299     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
300     } else {
301     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
302     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
303     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
304     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
305     }
306     }
307     }
308     totalNECount+=NE1*NE2;
309     faceNECount+=NE1*NE2;
310    
311     /* ** elements on boundary 002 (x1=1): */
312    
313     #pragma omp parallel for private(i1,i2,k,node0)
314     for (i2=0;i2<NE2;i2++) {
315     for (i1=0;i1<NE1;i1++) {
316     k=i1+NE1*i2+faceNECount;
317     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
318    
319     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
320     out->FaceElements->Tag[k]=2;
321     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
322    
323     if (useElementsOnFace) {
324     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
325     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
326     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
327     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
328     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
329     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
330     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
331     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
332     } else {
333     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
334     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
335     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
336     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
337     }
338     }
339     }
340     totalNECount+=NE1*NE2;
341     faceNECount+=NE1*NE2;
342     }
343     if (!periodic[1]) {
344     /* ** elements on boundary 010 (x2=0): */
345    
346     #pragma omp parallel for private(i0,i2,k,node0)
347     for (i2=0;i2<NE2;i2++) {
348     for (i0=0;i0<NE0;i0++) {
349     k=i0+NE0*i2+faceNECount;
350     node0=i0+N0*N1*i2;
351    
352     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
353     out->FaceElements->Tag[k]=10;
354     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
355    
356     if (useElementsOnFace) {
357     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
358     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
359     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
360     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
361     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
362     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
363     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
364     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
365     } else {
366     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
367     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
368     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
369     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
370     }
371     }
372     }
373     totalNECount+=NE0*NE2;
374     faceNECount+=NE0*NE2;
375    
376     /* ** elements on boundary 020 (x2=1): */
377    
378     #pragma omp parallel for private(i0,i2,k,node0)
379     for (i2=0;i2<NE2;i2++) {
380     for (i0=0;i0<NE0;i0++) {
381     k=i0+NE0*i2+faceNECount;
382     node0=i0+(NE1-1)*N0+N0*N1*i2;
383    
384     out->FaceElements->Tag[k]=20;
385     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
386     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
387    
388     if (useElementsOnFace) {
389     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
390     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
391     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
392     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
393     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
394     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
395     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
396     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
397     } else {
398     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
399     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
400     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
401     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
402     }
403    
404     }
405     }
406     totalNECount+=NE0*NE2;
407     faceNECount+=NE0*NE2;
408     }
409 jgs 123 out->FaceElements->minColor=0;
410     out->FaceElements->maxColor=23;
411 jgs 82
412     /* face elements done: */
413    
414     /* condense the nodes: */
415    
416     Finley_Mesh_resolveNodeIds(out);
417    
418     /* prepare mesh for further calculatuions:*/
419     Finley_Mesh_prepare(out) ;
420    
421 jgs 149 #ifdef Finley_TRACE
422 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
423 jgs 149 #endif
424 jgs 82
425 jgs 150 if (! Finley_noError()) {
426 jgs 82 Finley_Mesh_dealloc(out);
427     return NULL;
428     }
429     return out;
430     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26