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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 21953 byte(s)
Initial revision

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26