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

Annotation of /temp_trunk_copy/finley/src/Mesh_hex8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 471 - (hide annotations)
Fri Jan 27 01:33:02 2006 UTC (13 years, 9 months ago) by jgs
Original Path: trunk/finley/src/Mesh_hex8.c
File MIME type: text/plain
File size: 15104 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26