/[escript]/branches/arrayview_from_1695_trunk/finley/src/Mesh_hex20.c
ViewVC logotype

Contents of /branches/arrayview_from_1695_trunk/finley/src/Mesh_hex20.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 471 - (show annotations)
Fri Jan 27 01:33:02 2006 UTC (13 years, 3 months ago) by jgs
Original Path: trunk/finley/src/Mesh_hex20.c
File MIME type: text/plain
File size: 23143 byte(s)
reorganise finley src tree to remove inc dir and src/finley directory

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,M0,M1,M2;
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 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 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 if (! Finley_noError()) return NULL;
106
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 if (! Finley_noError()) {
117 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 if (! Finley_noError()) {
128 Finley_Mesh_dealloc(out);
129 return NULL;
130 }
131
132 #pragma omp parallel for private(i0,i1,i2,k)
133 for (i2=0;i2<N2;i2++) {
134 for (i1=0;i1<N1;i1++) {
135 for (i0=0;i0<N0;i0++) {
136 k=M0*i0+M1*i1+M2*i2;
137 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 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
141 out->Nodes->Tag[k]=0;
142 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
143 }
144 }
145 }
146
147
148 /* tags for the faces: */
149 if (!periodic[2]) {
150 for (i1=0;i1<N1;i1++) {
151 for (i0=0;i0<N0;i0++) {
152 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
153 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
154 }
155 }
156 }
157 if (!periodic[1]) {
158 for (i2=0;i2<N2;i2++) {
159 for (i0=0;i0<N0;i0++) {
160 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
161 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
162 }
163 }
164 }
165 if (!periodic[0]) {
166 for (i2=0;i2<N2;i2++) {
167 for (i1=0;i1<N1;i1++) {
168 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
169 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
170 }
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 out->Elements->minColor=0;
211 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
212
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 out->FaceElements->minColor=0;
540 out->FaceElements->maxColor=24;
541
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 #ifdef Finley_TRACE
552 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
553 #endif
554
555 if (! Finley_noError()) {
556 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