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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex8.c
File MIME type: text/plain
File size: 14479 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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 jgs 123 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
25     dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;
26     index_t node0;
27 jgs 82 Finley_Mesh* out;
28     char name[50];
29     double time0=Finley_timer();
30     NE0=MAX(1,numElements[0]);
31     NE1=MAX(1,numElements[1]);
32     NE2=MAX(1,numElements[2]);
33     N0=NE0+1;
34     N1=NE1+1;
35     N2=NE2+1;
36    
37     NFaceElements=0;
38     if (!periodic[0]) {
39     NDOF0=N0;
40     NFaceElements+=2*NE1*NE2;
41     } else {
42     NDOF0=N0-1;
43     }
44     if (!periodic[1]) {
45     NDOF1=N1;
46     NFaceElements+=2*NE0*NE2;
47     } else {
48     NDOF1=N1-1;
49     }
50     if (!periodic[2]) {
51     NDOF2=N2;
52     NFaceElements+=2*NE0*NE1;
53     } else {
54     NDOF2=N2-1;
55     }
56    
57     /* allocate mesh: */
58    
59     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
60     out=Finley_Mesh_alloc(name,3,order);
61     if (Finley_ErrorCode!=NO_ERROR) return NULL;
62    
63     out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
64     if (useElementsOnFace) {
65     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
66     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
67     } else {
68     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
69     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
70     }
71     out->Points=Finley_ElementFile_alloc(Point1,out->order);
72     if (Finley_ErrorCode!=NO_ERROR) {
73     Finley_Mesh_dealloc(out);
74     return NULL;
75     }
76    
77    
78     /* allocate tables: */
79    
80     Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
81     Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
82     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
83     if (Finley_ErrorCode!=NO_ERROR) {
84     Finley_Mesh_dealloc(out);
85     return NULL;
86     }
87    
88     /* set nodes: */
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     /* 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=i0+i1*N0+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,8)]=node0;
144     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
145     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
146     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
147     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
148     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
149     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
150     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
151    
152     }
153     }
154     }
155 jgs 123 out->Elements->minColor=0;
156     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
157 jgs 82
158     /* face elements: */
159    
160     if (useElementsOnFace) {
161     NUMNODES=8;
162     } else {
163     NUMNODES=4;
164     }
165     totalNECount=NE0*NE1*NE2;
166     faceNECount=0;
167    
168     /* these are the quadrilateral elements on boundary 1 (x3=0): */
169     if (!periodic[2]) {
170     /* ** elements on boundary 100 (x3=0): */
171    
172     #pragma omp parallel for private(i0,i1,k,node0)
173     for (i1=0;i1<NE1;i1++) {
174     for (i0=0;i0<NE0;i0++) {
175     k=i0+NE0*i1+faceNECount;
176     node0=i0+i1*N0;
177    
178     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
179     out->FaceElements->Tag[k]=100;
180     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
181    
182     if (useElementsOnFace) {
183     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
184     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
185     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
186     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
187     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
188     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
189     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
190     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
191     } else {
192     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
193     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
194     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
195     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
196     }
197     }
198     }
199     totalNECount+=NE1*NE0;
200     faceNECount+=NE1*NE0;
201    
202     /* ** elements on boundary 200 (x3=1) */
203    
204     #pragma omp parallel for private(i0,i1,k,node0)
205     for (i1=0;i1<NE1;i1++) {
206     for (i0=0;i0<NE0;i0++) {
207     k=i0+NE0*i1+faceNECount;
208     node0=i0+i1*N0+N0*N1*(NE2-1);
209    
210     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
211     out->FaceElements->Tag[k]=200;
212     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
213    
214     if (useElementsOnFace) {
215     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
216     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
217     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
218     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
219     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
220     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
221     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
222     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
223     } else {
224     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
225     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
226     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
227     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
228     }
229    
230    
231     }
232     }
233     totalNECount+=NE1*NE0;
234     faceNECount+=NE1*NE0;
235     }
236     if (!periodic[0]) {
237     /* ** elements on boundary 001 (x1=0): */
238    
239     #pragma omp parallel for private(i1,i2,k,node0)
240     for (i2=0;i2<NE2;i2++) {
241     for (i1=0;i1<NE1;i1++) {
242     k=i1+NE1*i2+faceNECount;
243     node0=i1*N0+N0*N1*i2;
244    
245     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
246     out->FaceElements->Tag[k]=1;
247     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
248    
249     if (useElementsOnFace) {
250     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
251     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
252     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
253     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
254     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
255     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
256     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
257     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
258     } else {
259     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
260     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
261     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
262     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
263     }
264     }
265     }
266     totalNECount+=NE1*NE2;
267     faceNECount+=NE1*NE2;
268    
269     /* ** elements on boundary 002 (x1=1): */
270    
271     #pragma omp parallel for private(i1,i2,k,node0)
272     for (i2=0;i2<NE2;i2++) {
273     for (i1=0;i1<NE1;i1++) {
274     k=i1+NE1*i2+faceNECount;
275     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
276    
277     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
278     out->FaceElements->Tag[k]=2;
279     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
280    
281     if (useElementsOnFace) {
282     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
283     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
284     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
285     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
286     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
287     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
288     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
289     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
290     } else {
291     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
292     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
293     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
294     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
295     }
296     }
297     }
298     totalNECount+=NE1*NE2;
299     faceNECount+=NE1*NE2;
300     }
301     if (!periodic[1]) {
302     /* ** elements on boundary 010 (x2=0): */
303    
304     #pragma omp parallel for private(i0,i2,k,node0)
305     for (i2=0;i2<NE2;i2++) {
306     for (i0=0;i0<NE0;i0++) {
307     k=i0+NE0*i2+faceNECount;
308     node0=i0+N0*N1*i2;
309    
310     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
311     out->FaceElements->Tag[k]=10;
312     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
313    
314     if (useElementsOnFace) {
315     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
316     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
317     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
318     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
319     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
320     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
321     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
322     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
323     } else {
324     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
325     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
326     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
327     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
328     }
329     }
330     }
331     totalNECount+=NE0*NE2;
332     faceNECount+=NE0*NE2;
333    
334     /* ** elements on boundary 020 (x2=1): */
335    
336     #pragma omp parallel for private(i0,i2,k,node0)
337     for (i2=0;i2<NE2;i2++) {
338     for (i0=0;i0<NE0;i0++) {
339     k=i0+NE0*i2+faceNECount;
340     node0=i0+(NE1-1)*N0+N0*N1*i2;
341    
342     out->FaceElements->Tag[k]=20;
343     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
344     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
345    
346     if (useElementsOnFace) {
347     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
348     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
349     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
350     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
351     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
352     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
353     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
354     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
355     } else {
356     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
357     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
358     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
359     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
360     }
361    
362     }
363     }
364     totalNECount+=NE0*NE2;
365     faceNECount+=NE0*NE2;
366     }
367 jgs 123 out->FaceElements->minColor=0;
368     out->FaceElements->maxColor=23;
369 jgs 82
370     /* face elements done: */
371    
372     /* condense the nodes: */
373    
374     Finley_Mesh_resolveNodeIds(out);
375    
376     /* prepare mesh for further calculatuions:*/
377     Finley_Mesh_prepare(out) ;
378    
379 jgs 149 #ifdef Finley_TRACE
380 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
381 jgs 149 #endif
382 jgs 82
383     if (Finley_ErrorCode!=NO_ERROR) {
384     Finley_Mesh_dealloc(out);
385     return NULL;
386     }
387     return out;
388     }
389    
390     /*
391     * $Log$
392 jgs 149 * Revision 1.3 2005/09/01 03:31:35 jgs
393     * Merge of development branch dev-02 back to main trunk on 2005-09-01
394     *
395     * Revision 1.2.2.1 2005/08/24 02:02:18 gross
396     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
397     *
398 jgs 123 * Revision 1.2 2005/07/08 04:07:52 jgs
399     * Merge of development branch back to main trunk on 2005-07-08
400 jgs 82 *
401 jgs 123 * Revision 1.1.1.1.2.1 2005/06/29 02:34:52 gross
402     * some changes towards 64 integers in finley
403     *
404     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
405     * initial import of project esys2
406     *
407 jgs 82 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
408     * Initial version of eys using boost-python.
409     *
410     *
411     */
412    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26