/[escript]/trunk/finley/src/Mesh_hex8.c
ViewVC logotype

Contents of /trunk/finley/src/Mesh_hex8.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: 15427 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 first order elements (Hex8) 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 /* Author: gross@access.edu.au */
24 /* Version: $Id$ */
25
26 /**************************************************************/
27
28 #include "RectangularMesh.h"
29
30 /**************************************************************/
31
32 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {
33 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
34 index_t node0;
35 Finley_Mesh* out;
36 char name[50];
37 double time0=Finley_timer();
38 NE0=MAX(1,numElements[0]);
39 NE1=MAX(1,numElements[1]);
40 NE2=MAX(1,numElements[2]);
41 N0=NE0+1;
42 N1=NE1+1;
43 N2=NE2+1;
44
45 if (N0<=MIN(N1,N2)) {
46 if (N1 <= N2) {
47 M0=1;
48 M1=N0;
49 M2=N0*N1;
50 } else {
51 M0=1;
52 M2=N0;
53 M1=N0*N2;
54 }
55 } else if (N1<=MIN(N2,N0)) {
56 if (N2 <= N0) {
57 M1=1;
58 M2=N1;
59 M0=N2*N1;
60 } else {
61 M1=1;
62 M0=N1;
63 M2=N1*N0;
64 }
65 } else {
66 if (N0 <= N1) {
67 M2=1;
68 M0=N2;
69 M1=N2*N0;
70 } else {
71 M2=1;
72 M1=N2;
73 M0=N1*N2;
74 }
75 }
76
77
78 NFaceElements=0;
79 if (!periodic[0]) {
80 NDOF0=N0;
81 NFaceElements+=2*NE1*NE2;
82 } else {
83 NDOF0=N0-1;
84 }
85 if (!periodic[1]) {
86 NDOF1=N1;
87 NFaceElements+=2*NE0*NE2;
88 } else {
89 NDOF1=N1-1;
90 }
91 if (!periodic[2]) {
92 NDOF2=N2;
93 NFaceElements+=2*NE0*NE1;
94 } else {
95 NDOF2=N2-1;
96 }
97
98 /* allocate mesh: */
99
100 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
101 /* TEMPFIX */
102 #ifndef PASO_MPI
103 out=Finley_Mesh_alloc(name,3,order);
104 #else
105 out=Finley_Mesh_alloc(name,3,order,NULL);
106 #endif
107 if (! Finley_noError()) return NULL;
108
109 #ifndef PASO_MPI
110 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
111 if (useElementsOnFace) {
112 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
113 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
114 } else {
115 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
116 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
117 }
118 out->Points=Finley_ElementFile_alloc(Point1,out->order);
119 #else
120 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,NULL);
121 if (useElementsOnFace) {
122 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,NULL);
123 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,NULL);
124 } else {
125 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,NULL);
126 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,NULL);
127 }
128 out->Points=Finley_ElementFile_alloc(Point1,out->order,NULL);
129 #endif
130 if (! Finley_noError()) {
131 Finley_Mesh_dealloc(out);
132 return NULL;
133 }
134
135
136 /* allocate tables: */
137 #ifndef PASO_MPI
138 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
139 #else
140 /* TODO */
141 #endif
142 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
143 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
144 if (! Finley_noError()) {
145 Finley_Mesh_dealloc(out);
146 return NULL;
147 }
148
149 /* set nodes: */
150
151 #pragma omp parallel for private(i0,i1,i2,k)
152 for (i2=0;i2<N2;i2++) {
153 for (i1=0;i1<N1;i1++) {
154 for (i0=0;i0<N0;i0++) {
155 k=M0*i0+M1*i1+M2*i2;
156 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
157 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
158 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
159 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
160 out->Nodes->Tag[k]=0;
161 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
162 }
163 }
164 }
165 /* tags for the faces: */
166 /* tags for the faces: */
167 if (!periodic[2]) {
168 for (i1=0;i1<N1;i1++) {
169 for (i0=0;i0<N0;i0++) {
170 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
171 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
172 }
173 }
174 }
175 if (!periodic[1]) {
176 for (i2=0;i2<N2;i2++) {
177 for (i0=0;i0<N0;i0++) {
178 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
179 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
180 }
181 }
182 }
183 if (!periodic[0]) {
184 for (i2=0;i2<N2;i2++) {
185 for (i1=0;i1<N1;i1++) {
186 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
187 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
188 }
189 }
190 }
191
192 /* set the elements: */
193
194 #pragma omp parallel for private(i0,i1,i2,k,node0)
195 for (i2=0;i2<NE2;i2++) {
196 for (i1=0;i1<NE1;i1++) {
197 for (i0=0;i0<NE0;i0++) {
198 k=i0+NE0*i1+NE0*NE1*i2;
199 node0=i0+i1*N0+N0*N1*i2;
200
201 out->Elements->Id[k]=k;
202 out->Elements->Tag[k]=0;
203 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
204
205 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
206 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
207 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
208 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
209 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
210 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
211 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
212 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
213
214 }
215 }
216 }
217 out->Elements->minColor=0;
218 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
219
220 /* face elements: */
221
222 if (useElementsOnFace) {
223 NUMNODES=8;
224 } else {
225 NUMNODES=4;
226 }
227 totalNECount=NE0*NE1*NE2;
228 faceNECount=0;
229
230 /* these are the quadrilateral elements on boundary 1 (x3=0): */
231 if (!periodic[2]) {
232 /* ** elements on boundary 100 (x3=0): */
233
234 #pragma omp parallel for private(i0,i1,k,node0)
235 for (i1=0;i1<NE1;i1++) {
236 for (i0=0;i0<NE0;i0++) {
237 k=i0+NE0*i1+faceNECount;
238 node0=i0+i1*N0;
239
240 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
241 out->FaceElements->Tag[k]=100;
242 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
243
244 if (useElementsOnFace) {
245 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
246 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
247 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
248 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
249 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
250 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
251 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
252 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
253 } else {
254 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
255 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
256 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
257 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
258 }
259 }
260 }
261 totalNECount+=NE1*NE0;
262 faceNECount+=NE1*NE0;
263
264 /* ** elements on boundary 200 (x3=1) */
265
266 #pragma omp parallel for private(i0,i1,k,node0)
267 for (i1=0;i1<NE1;i1++) {
268 for (i0=0;i0<NE0;i0++) {
269 k=i0+NE0*i1+faceNECount;
270 node0=i0+i1*N0+N0*N1*(NE2-1);
271
272 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
273 out->FaceElements->Tag[k]=200;
274 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
275
276 if (useElementsOnFace) {
277 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
278 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
279 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
280 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
281 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
282 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
283 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
284 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
285 } else {
286 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
287 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
288 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
289 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
290 }
291
292
293 }
294 }
295 totalNECount+=NE1*NE0;
296 faceNECount+=NE1*NE0;
297 }
298 if (!periodic[0]) {
299 /* ** elements on boundary 001 (x1=0): */
300
301 #pragma omp parallel for private(i1,i2,k,node0)
302 for (i2=0;i2<NE2;i2++) {
303 for (i1=0;i1<NE1;i1++) {
304 k=i1+NE1*i2+faceNECount;
305 node0=i1*N0+N0*N1*i2;
306
307 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
308 out->FaceElements->Tag[k]=1;
309 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
310
311 if (useElementsOnFace) {
312 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
313 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
314 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
315 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
316 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
317 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
318 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
319 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
320 } else {
321 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
322 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
323 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
324 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
325 }
326 }
327 }
328 totalNECount+=NE1*NE2;
329 faceNECount+=NE1*NE2;
330
331 /* ** elements on boundary 002 (x1=1): */
332
333 #pragma omp parallel for private(i1,i2,k,node0)
334 for (i2=0;i2<NE2;i2++) {
335 for (i1=0;i1<NE1;i1++) {
336 k=i1+NE1*i2+faceNECount;
337 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
338
339 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
340 out->FaceElements->Tag[k]=2;
341 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
342
343 if (useElementsOnFace) {
344 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
345 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
346 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
347 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
348 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
349 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
350 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
351 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
352 } else {
353 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
354 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
355 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
356 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
357 }
358 }
359 }
360 totalNECount+=NE1*NE2;
361 faceNECount+=NE1*NE2;
362 }
363 if (!periodic[1]) {
364 /* ** elements on boundary 010 (x2=0): */
365
366 #pragma omp parallel for private(i0,i2,k,node0)
367 for (i2=0;i2<NE2;i2++) {
368 for (i0=0;i0<NE0;i0++) {
369 k=i0+NE0*i2+faceNECount;
370 node0=i0+N0*N1*i2;
371
372 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
373 out->FaceElements->Tag[k]=10;
374 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
375
376 if (useElementsOnFace) {
377 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
378 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
379 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
380 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
381 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
382 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
383 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
384 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
385 } else {
386 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
387 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
388 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
389 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
390 }
391 }
392 }
393 totalNECount+=NE0*NE2;
394 faceNECount+=NE0*NE2;
395
396 /* ** elements on boundary 020 (x2=1): */
397
398 #pragma omp parallel for private(i0,i2,k,node0)
399 for (i2=0;i2<NE2;i2++) {
400 for (i0=0;i0<NE0;i0++) {
401 k=i0+NE0*i2+faceNECount;
402 node0=i0+(NE1-1)*N0+N0*N1*i2;
403
404 out->FaceElements->Tag[k]=20;
405 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
406 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
407
408 if (useElementsOnFace) {
409 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
410 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
411 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
412 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
413 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
414 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
415 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
416 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
417 } else {
418 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
419 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
420 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
421 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
422 }
423
424 }
425 }
426 totalNECount+=NE0*NE2;
427 faceNECount+=NE0*NE2;
428 }
429 out->FaceElements->minColor=0;
430 out->FaceElements->maxColor=23;
431
432 /* face elements done: */
433
434 /* condense the nodes: */
435
436 Finley_Mesh_resolveNodeIds(out);
437
438 /* prepare mesh for further calculatuions:*/
439 Finley_Mesh_prepare(out) ;
440
441 #ifdef Finley_TRACE
442 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
443 #endif
444
445 if (! Finley_noError()) {
446 Finley_Mesh_dealloc(out);
447 return NULL;
448 }
449 return out;
450 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26