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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 616 - (hide annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 5 months ago) by elspeth
File MIME type: text/plain
File size: 22821 byte(s)
Copyright added to more source files.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26