/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex20.c
File MIME type: text/plain
File size: 22259 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26