/[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 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 6 months ago) by elspeth
File MIME type: text/plain
File size: 22821 byte(s)
Copyright added to more source files.

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 out=Finley_Mesh_alloc(name,3,order);
103 if (! Finley_noError()) return NULL;
104
105 out->Elements=Finley_ElementFile_alloc(Hex20,out->order);
106 if (useElementsOnFace) {
107 out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order);
108 out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order);
109 } else {
110 out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order);
111 out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order);
112 }
113 out->Points=Finley_ElementFile_alloc(Point1,out->order);
114 if (! Finley_noError()) {
115 Finley_Mesh_dealloc(out);
116 return NULL;
117 }
118
119
120 /* allocate tables: */
121
122 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
123 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
124 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
125 if (! Finley_noError()) {
126 Finley_Mesh_dealloc(out);
127 return NULL;
128 }
129
130 #pragma omp parallel for private(i0,i1,i2,k)
131 for (i2=0;i2<N2;i2++) {
132 for (i1=0;i1<N1;i1++) {
133 for (i0=0;i0<N0;i0++) {
134 k=M0*i0+M1*i1+M2*i2;
135 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
136 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
137 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
138 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
139 out->Nodes->Tag[k]=0;
140 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
141 }
142 }
143 }
144
145
146 /* tags for the faces: */
147 if (!periodic[2]) {
148 for (i1=0;i1<N1;i1++) {
149 for (i0=0;i0<N0;i0++) {
150 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
151 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
152 }
153 }
154 }
155 if (!periodic[1]) {
156 for (i2=0;i2<N2;i2++) {
157 for (i0=0;i0<N0;i0++) {
158 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
159 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
160 }
161 }
162 }
163 if (!periodic[0]) {
164 for (i2=0;i2<N2;i2++) {
165 for (i1=0;i1<N1;i1++) {
166 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
167 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
168 }
169 }
170 }
171
172 /* set the elements: */
173
174 #pragma omp parallel for private(i0,i1,i2,k,node0)
175 for (i2=0;i2<NE2;i2++) {
176 for (i1=0;i1<NE1;i1++) {
177 for (i0=0;i0<NE0;i0++) {
178 k=i0+NE0*i1+NE0*NE1*i2;
179 node0=2*i0+2*i1*N0+2*N0*N1*i2;
180
181 out->Elements->Id[k]=k;
182 out->Elements->Tag[k]=0;
183 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
184
185 out->Elements->Nodes[INDEX2(0,k,20)]=node0;
186 out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
187 out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;
188 out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;
189 out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;
190 out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;
191 out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;
192 out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;
193 out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
194 out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;
195 out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;
196 out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;
197 out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;
198 out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;
199 out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;
200 out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;
201 out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;
202 out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;
203 out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;
204 out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;
205 }
206 }
207 }
208 out->Elements->minColor=0;
209 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
210
211 /* face elements: */
212
213 if (useElementsOnFace) {
214 NUMNODES=20;
215 } else {
216 NUMNODES=8;
217 }
218 totalNECount=NE0*NE1*NE2;
219 faceNECount=0;
220
221 /* these are the quadrilateral elements on boundary 1 (x3=0): */
222
223 if (!periodic[2]) {
224 /* ** elements on boundary 100 (x3=0): */
225 #pragma omp parallel for private(i0,i1,k,node0)
226 for (i1=0;i1<NE1;i1++) {
227 for (i0=0;i0<NE0;i0++) {
228 k=i0+NE0*i1+faceNECount;
229 node0=2*i0+2*i1*N0;
230
231 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
232 out->FaceElements->Tag[k]=100;
233 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
234
235 if (useElementsOnFace) {
236 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
237 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
238 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
239 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
240 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;
241 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;
242 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
243 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;
244 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;
245 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;
246 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;
247 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;
248 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
249 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;
250 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
251 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;
252 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;
253 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
254 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;
255 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;
256 } else {
257 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
258 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
259 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;
260 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;
261 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
262 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;
263 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;
264 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
265 }
266 }
267 }
268 totalNECount+=NE1*NE0;
269 faceNECount+=NE1*NE0;
270
271 /* ** elements on boundary 200 (x3=1) */
272 #pragma omp parallel for private(i0,i1,k,node0)
273 for (i1=0;i1<NE1;i1++) {
274 for (i0=0;i0<NE0;i0++) {
275 k=i0+NE0*i1+faceNECount;
276 node0=2*i0+2*i1*N0+2*N0*N1*(NE2-1);
277
278 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
279 out->FaceElements->Tag[k]=200;
280 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
281
282 if (useElementsOnFace) {
283 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
284 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
285 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
286 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
287
288 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
289 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;
290 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;
291 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;
292
293 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;
294 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;
295 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
296 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;
297
298 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;
299 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;
300 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;
301 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;
302
303 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;
304 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;
305 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;
306 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;
307
308 } else {
309 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
310 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;
311 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
312 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;
313 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;
314 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;
315 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
316 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;
317 }
318 }
319 }
320 totalNECount+=NE1*NE0;
321 faceNECount+=NE1*NE0;
322 }
323 if (!periodic[0]) {
324 /* ** elements on boundary 001 (x1=0): */
325
326 #pragma omp parallel for private(i1,i2,k,node0)
327 for (i2=0;i2<NE2;i2++) {
328 for (i1=0;i1<NE1;i1++) {
329 k=i1+NE1*i2+faceNECount;
330 node0=2*i1*N0+2*N0*N1*i2;
331
332 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
333 out->FaceElements->Tag[k]=1;
334 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
335
336 if (useElementsOnFace) {
337 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
338 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
339 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
340 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
341
342 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;
343 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;
344 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
345 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;
346
347 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;
348 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;
349 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;
350 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;
351
352 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
353 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;
354 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
355 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;
356
357 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;
358 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;
359 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;
360 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;
361
362 } else {
363 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
364 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;
365 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;
366 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;
367 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
368 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;
369 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;
370 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
371 }
372 }
373 }
374 totalNECount+=NE1*NE2;
375 faceNECount+=NE1*NE2;
376
377 /* ** elements on boundary 002 (x1=1): */
378
379 #pragma omp parallel for private(i1,i2,k,node0)
380 for (i2=0;i2<NE2;i2++) {
381 for (i1=0;i1<NE1;i1++) {
382 k=i1+NE1*i2+faceNECount;
383 node0=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;
384
385 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
386 out->FaceElements->Tag[k]=2;
387 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
388
389 if (useElementsOnFace) {
390 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
391 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
392 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
393 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
394
395 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
396 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;
397 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;
398 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;
399
400 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;
401 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;
402 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;
403 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;
404
405 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;
406 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;
407 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;
408 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;
409
410 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;
411 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;
412 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;
413 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;
414
415 } else {
416 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
417 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;
418 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
419 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;
420 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;
421 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;
422 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;
423 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;
424 }
425
426 }
427 }
428 totalNECount+=NE1*NE2;
429 faceNECount+=NE1*NE2;
430 }
431 if (!periodic[1]) {
432 /* ** elements on boundary 010 (x2=0): */
433
434 #pragma omp parallel for private(i0,i2,k,node0)
435 for (i2=0;i2<NE2;i2++) {
436 for (i0=0;i0<NE0;i0++) {
437 k=i0+NE0*i2+faceNECount;
438 node0=2*i0+2*N0*N1*i2;
439
440 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
441 out->FaceElements->Tag[k]=10;
442 out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
443
444 if (useElementsOnFace) {
445 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
446 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
447 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
448 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
449
450 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;
451 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;
452 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;
453 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;
454
455 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;
456 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;
457 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;
458 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;
459
460 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
461 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;
462 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;
463 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;
464
465 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;
466 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;
467 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
468 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;
469
470 } else {
471 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
472 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;
473 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;
474 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;
475 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
476 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;
477 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;
478 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;
479 }
480 }
481 }
482 totalNECount+=NE0*NE2;
483 faceNECount+=NE0*NE2;
484
485 /* ** elements on boundary 020 (x2=1): */
486
487 #pragma omp parallel for private(i0,i2,k,node0)
488 for (i2=0;i2<NE2;i2++) {
489 for (i0=0;i0<NE0;i0++) {
490 k=i0+NE0*i2+faceNECount;
491 node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;
492
493 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
494 out->FaceElements->Tag[k]=20;
495 out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
496
497 if (useElementsOnFace) {
498 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
499 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
500 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
501 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
502
503 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
504 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;
505 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;
506 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;
507
508 out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;
509 out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
510 out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;
511 out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;
512
513 out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;
514 out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;
515 out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;
516 out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;
517
518 out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;
519 out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;
520 out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;
521 out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;
522 } else {
523 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
524 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;
525 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;
526 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;
527 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;
528 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;
529 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;
530 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;
531 }
532 }
533 }
534 totalNECount+=NE0*NE2;
535 faceNECount+=NE0*NE2;
536 }
537 out->FaceElements->minColor=0;
538 out->FaceElements->maxColor=24;
539
540 /* face elements done: */
541
542 /* condense the nodes: */
543
544 Finley_Mesh_resolveNodeIds(out);
545
546 /* prepare mesh for further calculatuions:*/
547 Finley_Mesh_prepare(out) ;
548
549 #ifdef Finley_TRACE
550 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
551 #endif
552
553 if (! Finley_noError()) {
554 Finley_Mesh_dealloc(out);
555 return NULL;
556 }
557 return out;
558 }
559

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26