/[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 730 - (hide annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 5 months ago) by bcumming
File MIME type: text/plain
File size: 15427 byte(s)


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 bcumming 730 /* TEMPFIX */
102     #ifndef PASO_MPI
103 jgs 82 out=Finley_Mesh_alloc(name,3,order);
104 bcumming 730 #else
105     out=Finley_Mesh_alloc(name,3,order,NULL);
106     #endif
107 jgs 150 if (! Finley_noError()) return NULL;
108 jgs 82
109 bcumming 730 #ifndef PASO_MPI
110 jgs 82 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
111     if (useElementsOnFace) {
112     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
113     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
114     } else {
115     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
116     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
117     }
118     out->Points=Finley_ElementFile_alloc(Point1,out->order);
119 bcumming 730 #else
120     out->Elements=Finley_ElementFile_alloc(Hex8,out->order,NULL);
121     if (useElementsOnFace) {
122     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,NULL);
123     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,NULL);
124     } else {
125     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,NULL);
126     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,NULL);
127     }
128     out->Points=Finley_ElementFile_alloc(Point1,out->order,NULL);
129     #endif
130 jgs 150 if (! Finley_noError()) {
131 jgs 82 Finley_Mesh_dealloc(out);
132     return NULL;
133     }
134    
135    
136     /* allocate tables: */
137 bcumming 730 #ifndef PASO_MPI
138 jgs 82 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
139 bcumming 730 #else
140     /* TODO */
141     #endif
142 jgs 82 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
143     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
144 jgs 150 if (! Finley_noError()) {
145 jgs 82 Finley_Mesh_dealloc(out);
146     return NULL;
147     }
148    
149     /* set nodes: */
150 jgs 153
151     #pragma omp parallel for private(i0,i1,i2,k)
152 jgs 82 for (i2=0;i2<N2;i2++) {
153     for (i1=0;i1<N1;i1++) {
154     for (i0=0;i0<N0;i0++) {
155 jgs 153 k=M0*i0+M1*i1+M2*i2;
156 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
157     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
158     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
159 jgs 153 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
160 jgs 82 out->Nodes->Tag[k]=0;
161 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
162 jgs 82 }
163     }
164     }
165     /* tags for the faces: */
166 jgs 153 /* tags for the faces: */
167 jgs 82 if (!periodic[2]) {
168 jgs 153 for (i1=0;i1<N1;i1++) {
169     for (i0=0;i0<N0;i0++) {
170     out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
171     out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
172     }
173     }
174 jgs 82 }
175     if (!periodic[1]) {
176     for (i2=0;i2<N2;i2++) {
177     for (i0=0;i0<N0;i0++) {
178 jgs 153 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
179     out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
180 jgs 82 }
181     }
182     }
183     if (!periodic[0]) {
184     for (i2=0;i2<N2;i2++) {
185     for (i1=0;i1<N1;i1++) {
186 jgs 153 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
187     out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
188 jgs 82 }
189     }
190     }
191 jgs 153
192 jgs 82 /* set the elements: */
193    
194     #pragma omp parallel for private(i0,i1,i2,k,node0)
195     for (i2=0;i2<NE2;i2++) {
196     for (i1=0;i1<NE1;i1++) {
197     for (i0=0;i0<NE0;i0++) {
198     k=i0+NE0*i1+NE0*NE1*i2;
199     node0=i0+i1*N0+N0*N1*i2;
200    
201     out->Elements->Id[k]=k;
202     out->Elements->Tag[k]=0;
203     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
204    
205     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
206     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
207     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
208     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
209     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
210     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
211     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
212     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
213    
214     }
215     }
216     }
217 jgs 123 out->Elements->minColor=0;
218     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
219 jgs 82
220     /* face elements: */
221    
222     if (useElementsOnFace) {
223     NUMNODES=8;
224     } else {
225     NUMNODES=4;
226     }
227     totalNECount=NE0*NE1*NE2;
228     faceNECount=0;
229    
230     /* these are the quadrilateral elements on boundary 1 (x3=0): */
231     if (!periodic[2]) {
232     /* ** elements on boundary 100 (x3=0): */
233    
234     #pragma omp parallel for private(i0,i1,k,node0)
235     for (i1=0;i1<NE1;i1++) {
236     for (i0=0;i0<NE0;i0++) {
237     k=i0+NE0*i1+faceNECount;
238     node0=i0+i1*N0;
239    
240     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
241     out->FaceElements->Tag[k]=100;
242     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
243    
244     if (useElementsOnFace) {
245     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
246     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
247     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
248     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
249     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
250     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
251     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
252     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
253     } else {
254     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
255     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
256     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
257     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
258     }
259     }
260     }
261     totalNECount+=NE1*NE0;
262     faceNECount+=NE1*NE0;
263    
264     /* ** elements on boundary 200 (x3=1) */
265    
266     #pragma omp parallel for private(i0,i1,k,node0)
267     for (i1=0;i1<NE1;i1++) {
268     for (i0=0;i0<NE0;i0++) {
269     k=i0+NE0*i1+faceNECount;
270     node0=i0+i1*N0+N0*N1*(NE2-1);
271    
272     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
273     out->FaceElements->Tag[k]=200;
274     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
275    
276     if (useElementsOnFace) {
277     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
278     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
279     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
280     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
281     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
282     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
283     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
284     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
285     } else {
286     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
287     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
288     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
289     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
290     }
291    
292    
293     }
294     }
295     totalNECount+=NE1*NE0;
296     faceNECount+=NE1*NE0;
297     }
298     if (!periodic[0]) {
299     /* ** elements on boundary 001 (x1=0): */
300    
301     #pragma omp parallel for private(i1,i2,k,node0)
302     for (i2=0;i2<NE2;i2++) {
303     for (i1=0;i1<NE1;i1++) {
304     k=i1+NE1*i2+faceNECount;
305     node0=i1*N0+N0*N1*i2;
306    
307     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
308     out->FaceElements->Tag[k]=1;
309     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
310    
311     if (useElementsOnFace) {
312     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
313     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
314     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
315     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
316     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
317     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
318     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
319     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
320     } else {
321     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
322     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
323     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
324     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
325     }
326     }
327     }
328     totalNECount+=NE1*NE2;
329     faceNECount+=NE1*NE2;
330    
331     /* ** elements on boundary 002 (x1=1): */
332    
333     #pragma omp parallel for private(i1,i2,k,node0)
334     for (i2=0;i2<NE2;i2++) {
335     for (i1=0;i1<NE1;i1++) {
336     k=i1+NE1*i2+faceNECount;
337     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
338    
339     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
340     out->FaceElements->Tag[k]=2;
341     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
342    
343     if (useElementsOnFace) {
344     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
345     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
346     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
347     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
348     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
349     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
350     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
351     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
352     } else {
353     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
354     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
355     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
356     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
357     }
358     }
359     }
360     totalNECount+=NE1*NE2;
361     faceNECount+=NE1*NE2;
362     }
363     if (!periodic[1]) {
364     /* ** elements on boundary 010 (x2=0): */
365    
366     #pragma omp parallel for private(i0,i2,k,node0)
367     for (i2=0;i2<NE2;i2++) {
368     for (i0=0;i0<NE0;i0++) {
369     k=i0+NE0*i2+faceNECount;
370     node0=i0+N0*N1*i2;
371    
372     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
373     out->FaceElements->Tag[k]=10;
374     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
375    
376     if (useElementsOnFace) {
377     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
378     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
379     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
380     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
381     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
382     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
383     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
384     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
385     } else {
386     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
387     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
388     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
389     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
390     }
391     }
392     }
393     totalNECount+=NE0*NE2;
394     faceNECount+=NE0*NE2;
395    
396     /* ** elements on boundary 020 (x2=1): */
397    
398     #pragma omp parallel for private(i0,i2,k,node0)
399     for (i2=0;i2<NE2;i2++) {
400     for (i0=0;i0<NE0;i0++) {
401     k=i0+NE0*i2+faceNECount;
402     node0=i0+(NE1-1)*N0+N0*N1*i2;
403    
404     out->FaceElements->Tag[k]=20;
405     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
406     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
407    
408     if (useElementsOnFace) {
409     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
410     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
411     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
412     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
413     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
414     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
415     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
416     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
417     } else {
418     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
419     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
420     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
421     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
422     }
423    
424     }
425     }
426     totalNECount+=NE0*NE2;
427     faceNECount+=NE0*NE2;
428     }
429 jgs 123 out->FaceElements->minColor=0;
430     out->FaceElements->maxColor=23;
431 jgs 82
432     /* face elements done: */
433    
434     /* condense the nodes: */
435    
436     Finley_Mesh_resolveNodeIds(out);
437    
438     /* prepare mesh for further calculatuions:*/
439     Finley_Mesh_prepare(out) ;
440    
441 jgs 149 #ifdef Finley_TRACE
442 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
443 jgs 149 #endif
444 jgs 82
445 jgs 150 if (! Finley_noError()) {
446 jgs 82 Finley_Mesh_dealloc(out);
447     return NULL;
448     }
449     return out;
450     }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26