/[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 730 - (show annotations)
Mon May 15 04:03:49 2006 UTC (13 years, 5 months ago) by bcumming
File MIME type: text/plain
File size: 22952 byte(s)


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26