/[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 150 - (show annotations)
Thu Sep 15 03:44:45 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: 23406 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26