/[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 201 - (show annotations)
Wed Nov 23 04:10:21 2005 UTC (13 years, 9 months ago) by jgs
Original Path: trunk/finley/src/finley/Mesh_hex8.c
File MIME type: text/plain
File size: 15104 byte(s)
copy finleyC and CPPAdapter to finley and finley/CPPAdapter to
facilitate scons builds

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26