/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex8.c
File MIME type: text/plain
File size: 15338 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26