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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26