/[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 616 - (show annotations)
Wed Mar 22 02:46:56 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/plain
File size: 14782 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 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 out=Finley_Mesh_alloc(name,3,order);
102 if (! Finley_noError()) return NULL;
103
104 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
105 if (useElementsOnFace) {
106 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
107 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
108 } else {
109 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
110 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
111 }
112 out->Points=Finley_ElementFile_alloc(Point1,out->order);
113 if (! Finley_noError()) {
114 Finley_Mesh_dealloc(out);
115 return NULL;
116 }
117
118
119 /* allocate tables: */
120
121 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
122 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
123 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
124 if (! Finley_noError()) {
125 Finley_Mesh_dealloc(out);
126 return NULL;
127 }
128
129 /* set nodes: */
130
131 #pragma omp parallel for private(i0,i1,i2,k)
132 for (i2=0;i2<N2;i2++) {
133 for (i1=0;i1<N1;i1++) {
134 for (i0=0;i0<N0;i0++) {
135 k=M0*i0+M1*i1+M2*i2;
136 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
137 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
138 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
139 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
140 out->Nodes->Tag[k]=0;
141 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
142 }
143 }
144 }
145 /* tags for the faces: */
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=i0+i1*N0+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,8)]=node0;
186 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
187 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
188 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
189 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
190 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
191 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
192 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
193
194 }
195 }
196 }
197 out->Elements->minColor=0;
198 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
199
200 /* face elements: */
201
202 if (useElementsOnFace) {
203 NUMNODES=8;
204 } else {
205 NUMNODES=4;
206 }
207 totalNECount=NE0*NE1*NE2;
208 faceNECount=0;
209
210 /* these are the quadrilateral elements on boundary 1 (x3=0): */
211 if (!periodic[2]) {
212 /* ** elements on boundary 100 (x3=0): */
213
214 #pragma omp parallel for private(i0,i1,k,node0)
215 for (i1=0;i1<NE1;i1++) {
216 for (i0=0;i0<NE0;i0++) {
217 k=i0+NE0*i1+faceNECount;
218 node0=i0+i1*N0;
219
220 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
221 out->FaceElements->Tag[k]=100;
222 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
223
224 if (useElementsOnFace) {
225 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
226 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
227 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
228 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
229 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
230 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
231 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
232 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
233 } else {
234 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
235 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
236 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
237 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
238 }
239 }
240 }
241 totalNECount+=NE1*NE0;
242 faceNECount+=NE1*NE0;
243
244 /* ** elements on boundary 200 (x3=1) */
245
246 #pragma omp parallel for private(i0,i1,k,node0)
247 for (i1=0;i1<NE1;i1++) {
248 for (i0=0;i0<NE0;i0++) {
249 k=i0+NE0*i1+faceNECount;
250 node0=i0+i1*N0+N0*N1*(NE2-1);
251
252 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
253 out->FaceElements->Tag[k]=200;
254 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
255
256 if (useElementsOnFace) {
257 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
258 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
259 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
260 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
261 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
262 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
263 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
264 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
265 } else {
266 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
267 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
268 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
269 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
270 }
271
272
273 }
274 }
275 totalNECount+=NE1*NE0;
276 faceNECount+=NE1*NE0;
277 }
278 if (!periodic[0]) {
279 /* ** elements on boundary 001 (x1=0): */
280
281 #pragma omp parallel for private(i1,i2,k,node0)
282 for (i2=0;i2<NE2;i2++) {
283 for (i1=0;i1<NE1;i1++) {
284 k=i1+NE1*i2+faceNECount;
285 node0=i1*N0+N0*N1*i2;
286
287 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
288 out->FaceElements->Tag[k]=1;
289 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
290
291 if (useElementsOnFace) {
292 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
293 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
294 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
295 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
296 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
297 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
298 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
299 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
300 } else {
301 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
302 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
303 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
304 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
305 }
306 }
307 }
308 totalNECount+=NE1*NE2;
309 faceNECount+=NE1*NE2;
310
311 /* ** elements on boundary 002 (x1=1): */
312
313 #pragma omp parallel for private(i1,i2,k,node0)
314 for (i2=0;i2<NE2;i2++) {
315 for (i1=0;i1<NE1;i1++) {
316 k=i1+NE1*i2+faceNECount;
317 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
318
319 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
320 out->FaceElements->Tag[k]=2;
321 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
322
323 if (useElementsOnFace) {
324 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
325 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
326 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
327 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
328 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
329 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
330 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
331 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
332 } else {
333 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
334 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
335 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
336 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
337 }
338 }
339 }
340 totalNECount+=NE1*NE2;
341 faceNECount+=NE1*NE2;
342 }
343 if (!periodic[1]) {
344 /* ** elements on boundary 010 (x2=0): */
345
346 #pragma omp parallel for private(i0,i2,k,node0)
347 for (i2=0;i2<NE2;i2++) {
348 for (i0=0;i0<NE0;i0++) {
349 k=i0+NE0*i2+faceNECount;
350 node0=i0+N0*N1*i2;
351
352 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
353 out->FaceElements->Tag[k]=10;
354 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
355
356 if (useElementsOnFace) {
357 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
358 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
359 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
360 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
361 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
362 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
363 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
364 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
365 } else {
366 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
367 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
368 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
369 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
370 }
371 }
372 }
373 totalNECount+=NE0*NE2;
374 faceNECount+=NE0*NE2;
375
376 /* ** elements on boundary 020 (x2=1): */
377
378 #pragma omp parallel for private(i0,i2,k,node0)
379 for (i2=0;i2<NE2;i2++) {
380 for (i0=0;i0<NE0;i0++) {
381 k=i0+NE0*i2+faceNECount;
382 node0=i0+(NE1-1)*N0+N0*N1*i2;
383
384 out->FaceElements->Tag[k]=20;
385 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
386 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
387
388 if (useElementsOnFace) {
389 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
390 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
391 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
392 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
393 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
394 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
395 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
396 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
397 } else {
398 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
399 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
400 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
401 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
402 }
403
404 }
405 }
406 totalNECount+=NE0*NE2;
407 faceNECount+=NE0*NE2;
408 }
409 out->FaceElements->minColor=0;
410 out->FaceElements->maxColor=23;
411
412 /* face elements done: */
413
414 /* condense the nodes: */
415
416 Finley_Mesh_resolveNodeIds(out);
417
418 /* prepare mesh for further calculatuions:*/
419 Finley_Mesh_prepare(out) ;
420
421 #ifdef Finley_TRACE
422 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
423 #endif
424
425 if (! Finley_noError()) {
426 Finley_Mesh_dealloc(out);
427 return NULL;
428 }
429 return out;
430 }

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26