/[escript]/branches/domexper/dudley/src/Mesh_hex20.c
ViewVC logotype

Annotation of /branches/domexper/dudley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 5 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 23406 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 second order elements (Hex20) 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    
26 jgs 150 /* Author: gross@access.edu.au */
27     /* Version: $Id$
28 jgs 82
29     /**************************************************************/
30    
31     #include "RectangularMesh.h"
32    
33     /**************************************************************/
34    
35 jgs 123 Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {
36     dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;
37     index_t node0;
38 jgs 82 Finley_Mesh* out;
39     char name[50];
40     double time0=Finley_timer();
41    
42     NE0=MAX(1,numElements[0]);
43     NE1=MAX(1,numElements[1]);
44     NE2=MAX(1,numElements[2]);
45     N0=2*NE0+1;
46     N1=2*NE1+1;
47     N2=2*NE2+1;
48    
49     NFaceElements=0;
50     if (!periodic[0]) {
51     NDOF0=N0;
52     NFaceElements+=2*NE1*NE2;
53     } else {
54     NDOF0=N0-1;
55     }
56     if (!periodic[1]) {
57     NDOF1=N1;
58     NFaceElements+=2*NE0*NE2;
59     } else {
60     NDOF1=N1-1;
61     }
62     if (!periodic[2]) {
63     NDOF2=N2;
64     NFaceElements+=2*NE0*NE1;
65     } else {
66     NDOF2=N2-1;
67     }
68    
69     /* allocate mesh: */
70    
71     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
72     out=Finley_Mesh_alloc(name,3,order);
73 jgs 150 if (! Finley_noError()) return NULL;
74 jgs 82
75     out->Elements=Finley_ElementFile_alloc(Hex20,out->order);
76     if (useElementsOnFace) {
77     out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order);
78     out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order);
79     } else {
80     out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order);
81     out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order);
82     }
83     out->Points=Finley_ElementFile_alloc(Point1,out->order);
84 jgs 150 if (! Finley_noError()) {
85 jgs 82 Finley_Mesh_dealloc(out);
86     return NULL;
87     }
88    
89    
90     /* allocate tables: */
91    
92     Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
93     Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
94     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
95 jgs 150 if (! Finley_noError()) {
96 jgs 82 Finley_Mesh_dealloc(out);
97     return NULL;
98     }
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    
115     /* tags for the faces: */
116     if (!periodic[2]) {
117     for (i1=0;i1<N1;i1++) {
118     for (i0=0;i0<N0;i0++) {
119     out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;
120     out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;
121     }
122     }
123     }
124     if (!periodic[1]) {
125     for (i2=0;i2<N2;i2++) {
126     for (i0=0;i0<N0;i0++) {
127     out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;
128     out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;
129     }
130     }
131     }
132     if (!periodic[0]) {
133     for (i2=0;i2<N2;i2++) {
134     for (i1=0;i1<N1;i1++) {
135     out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;
136     out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;
137     }
138     }
139     }
140    
141     /* set the elements: */
142    
143     #pragma omp parallel for private(i0,i1,i2,k,node0)
144     for (i2=0;i2<NE2;i2++) {
145     for (i1=0;i1<NE1;i1++) {
146     for (i0=0;i0<NE0;i0++) {
147     k=i0+NE0*i1+NE0*NE1*i2;
148     node0=2*i0+2*i1*N0+2*N0*N1*i2;
149    
150     out->Elements->Id[k]=k;
151     out->Elements->Tag[k]=0;
152     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
153    
154     out->Elements->Nodes[INDEX2(0,k,20)]=node0;
155     out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
156     out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;
157     out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;
158     out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;
159     out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;
160     out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;
161     out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;
162     out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
163     out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;
164     out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;
165     out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;
166     out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;
167     out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;
168     out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;
169     out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;
170     out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;
171     out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;
172     out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;
173     out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;
174     }
175     }
176     }
177 jgs 123 out->Elements->minColor=0;
178     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
179 jgs 82
180     /* face elements: */
181    
182     if (useElementsOnFace) {
183     NUMNODES=20;
184     } else {
185     NUMNODES=8;
186     }
187     totalNECount=NE0*NE1*NE2;
188     faceNECount=0;
189    
190     /* these are the quadrilateral elements on boundary 1 (x3=0): */
191    
192     if (!periodic[2]) {
193     /* ** elements on boundary 100 (x3=0): */
194     #pragma omp parallel for private(i0,i1,k,node0)
195     for (i1=0;i1<NE1;i1++) {
196     for (i0=0;i0<NE0;i0++) {
197     k=i0+NE0*i1+faceNECount;
198     node0=2*i0+2*i1*N0;
199    
200     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
201     out->FaceElements->Tag[k]=100;
202     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
203    
204     if (useElementsOnFace) {
205     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
206     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
207     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
208     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
209     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;
210     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;
211     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
212     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;
213     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;
214     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;
215     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;
216     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;
217     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
218     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;
219     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
220     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;
221     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;
222     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
223     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;
224     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;
225     } else {
226     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
227     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
228     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
229     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
230     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
231     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
232     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
233     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
234     }
235     }
236     }
237     totalNECount+=NE1*NE0;
238     faceNECount+=NE1*NE0;
239    
240     /* ** elements on boundary 200 (x3=1) */
241     #pragma omp parallel for private(i0,i1,k,node0)
242     for (i1=0;i1<NE1;i1++) {
243     for (i0=0;i0<NE0;i0++) {
244     k=i0+NE0*i1+faceNECount;
245     node0=2*i0+2*i1*N0+2*N0*N1*(NE2-1);
246    
247     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
248     out->FaceElements->Tag[k]=200;
249     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
250    
251     if (useElementsOnFace) {
252     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
253     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
254     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
255     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
256    
257     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
258     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;
259     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;
260     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;
261    
262     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;
263     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;
264     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
265     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;
266    
267     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
268     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;
269     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
270     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;
271    
272     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;
273     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;
274     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;
275     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;
276    
277     } else {
278     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
279     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
280     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
281     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
282     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;
283     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;
284     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
285     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;
286     }
287     }
288     }
289     totalNECount+=NE1*NE0;
290     faceNECount+=NE1*NE0;
291     }
292     if (!periodic[0]) {
293     /* ** elements on boundary 001 (x1=0): */
294    
295     #pragma omp parallel for private(i1,i2,k,node0)
296     for (i2=0;i2<NE2;i2++) {
297     for (i1=0;i1<NE1;i1++) {
298     k=i1+NE1*i2+faceNECount;
299     node0=2*i1*N0+2*N0*N1*i2;
300    
301     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
302     out->FaceElements->Tag[k]=1;
303     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
304    
305     if (useElementsOnFace) {
306     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
307     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
308     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
309     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
310    
311     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;
312     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;
313     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
314     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;
315    
316     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;
317     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;
318     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;
319     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;
320    
321     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
322     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;
323     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
324     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;
325    
326     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;
327     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;
328     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;
329     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;
330    
331     } else {
332     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
333     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
334     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
335     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
336     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
337     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;
338     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;
339     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
340     }
341     }
342     }
343     totalNECount+=NE1*NE2;
344     faceNECount+=NE1*NE2;
345    
346     /* ** elements on boundary 002 (x1=1): */
347    
348     #pragma omp parallel for private(i1,i2,k,node0)
349     for (i2=0;i2<NE2;i2++) {
350     for (i1=0;i1<NE1;i1++) {
351     k=i1+NE1*i2+faceNECount;
352     node0=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;
353    
354     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
355     out->FaceElements->Tag[k]=2;
356     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
357    
358     if (useElementsOnFace) {
359     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
360     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
361     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
362     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
363    
364     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
365     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;
366     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;
367     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;
368    
369     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;
370     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;
371     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;
372     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;
373    
374     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
375     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;
376     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
377     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;
378    
379     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;
380     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;
381     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;
382     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;
383    
384     } else {
385     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
386     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
387     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
388     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
389     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
390     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;
391     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;
392     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;
393     }
394    
395     }
396     }
397     totalNECount+=NE1*NE2;
398     faceNECount+=NE1*NE2;
399     }
400     if (!periodic[1]) {
401     /* ** elements on boundary 010 (x2=0): */
402    
403     #pragma omp parallel for private(i0,i2,k,node0)
404     for (i2=0;i2<NE2;i2++) {
405     for (i0=0;i0<NE0;i0++) {
406     k=i0+NE0*i2+faceNECount;
407     node0=2*i0+2*N0*N1*i2;
408    
409     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
410     out->FaceElements->Tag[k]=10;
411     out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
412    
413     if (useElementsOnFace) {
414     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
415     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
416     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
417     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
418    
419     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;
420     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;
421     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;
422     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;
423    
424     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;
425     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;
426     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;
427     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;
428    
429     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
430     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;
431     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;
432     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;
433    
434     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;
435     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;
436     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
437     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;
438    
439     } else {
440     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
441     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
442     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
443     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
444     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
445     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;
446     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;
447     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;
448     }
449     }
450     }
451     totalNECount+=NE0*NE2;
452     faceNECount+=NE0*NE2;
453    
454     /* ** elements on boundary 020 (x2=1): */
455    
456     #pragma omp parallel for private(i0,i2,k,node0)
457     for (i2=0;i2<NE2;i2++) {
458     for (i0=0;i0<NE0;i0++) {
459     k=i0+NE0*i2+faceNECount;
460     node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;
461    
462     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
463     out->FaceElements->Tag[k]=20;
464     out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
465    
466     if (useElementsOnFace) {
467     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
468     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
469     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
470     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
471    
472     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
473     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;
474     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;
475     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;
476    
477     out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;
478     out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
479     out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;
480     out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;
481    
482     out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
483     out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;
484     out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;
485     out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;
486    
487     out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;
488     out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;
489     out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;
490     out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;
491     } else {
492     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
493     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
494     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
495     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
496     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;
497     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
498     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;
499     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
500     }
501     }
502     }
503     totalNECount+=NE0*NE2;
504     faceNECount+=NE0*NE2;
505     }
506 jgs 123 out->FaceElements->minColor=0;
507     out->FaceElements->maxColor=24;
508 jgs 82
509     /* face elements done: */
510    
511     /* condense the nodes: */
512    
513     Finley_Mesh_resolveNodeIds(out);
514    
515     /* prepare mesh for further calculatuions:*/
516     Finley_Mesh_prepare(out) ;
517    
518 jgs 149 #ifdef Finley_TRACE
519 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
520 jgs 149 #endif
521 jgs 82
522 jgs 150 if (! Finley_noError()) {
523 jgs 82 Finley_Mesh_dealloc(out);
524     return NULL;
525     }
526     return out;
527     }
528    
529     /*
530 jgs 149 * Revision 1.3 2005/09/01 03:31:35 jgs
531     * Merge of development branch dev-02 back to main trunk on 2005-09-01
532     *
533 jgs 150 * Revision 1.2.2.2 2005/09/07 06:26:19 gross
534     * the solver from finley are put into the standalone package paso now
535     *
536 jgs 149 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
537     * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
538     *
539 jgs 123 * Revision 1.2 2005/07/08 04:07:52 jgs
540     * Merge of development branch back to main trunk on 2005-07-08
541 jgs 82 *
542 jgs 123 * Revision 1.1.1.1.2.1 2005/06/29 02:34:51 gross
543     * some changes towards 64 integers in finley
544     *
545     * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
546     * initial import of project esys2
547     *
548 jgs 82 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
549     * Initial version of eys using boost-python.
550     *
551     *
552     */
553    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26