/[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 153 - (hide annotations)
Tue Oct 25 01:51:20 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 23143 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-10-25

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26