/[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 82 - (show annotations)
Tue Oct 26 06:53:54 2004 UTC (15 years ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex8.c
File MIME type: text/plain
File size: 13882 byte(s)
Initial revision

1 /**************************************************************/
2
3 /* Finley: generates rectangular meshes */
4
5 /* Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
6 /* [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
7 /* integration scheme. */
8
9 /**************************************************************/
10
11 /* Copyrights by ACcESS Australia 2003/04 */
12 /* Author: gross@access.edu.au */
13 /* Version: $Id$ */
14
15 /**************************************************************/
16
17 #include "Common.h"
18 #include "Finley.h"
19 #include "Mesh.h"
20 #include "RectangularMesh.h"
21
22 /**************************************************************/
23
24 Finley_Mesh* Finley_RectangularMesh_Hex8(int* numElements,double* Length,int* periodic, int order,int useElementsOnFace) {
25 int N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,node0,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;
26 Finley_Mesh* out;
27 char name[50];
28 double time0=Finley_timer();
29 NE0=MAX(1,numElements[0]);
30 NE1=MAX(1,numElements[1]);
31 NE2=MAX(1,numElements[2]);
32 N0=NE0+1;
33 N1=NE1+1;
34 N2=NE2+1;
35
36 NFaceElements=0;
37 if (!periodic[0]) {
38 NDOF0=N0;
39 NFaceElements+=2*NE1*NE2;
40 } else {
41 NDOF0=N0-1;
42 }
43 if (!periodic[1]) {
44 NDOF1=N1;
45 NFaceElements+=2*NE0*NE2;
46 } else {
47 NDOF1=N1-1;
48 }
49 if (!periodic[2]) {
50 NDOF2=N2;
51 NFaceElements+=2*NE0*NE1;
52 } else {
53 NDOF2=N2-1;
54 }
55
56 /* allocate mesh: */
57
58 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
59 out=Finley_Mesh_alloc(name,3,order);
60 if (Finley_ErrorCode!=NO_ERROR) return NULL;
61
62 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
63 if (useElementsOnFace) {
64 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
65 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
66 } else {
67 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
68 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
69 }
70 out->Points=Finley_ElementFile_alloc(Point1,out->order);
71 if (Finley_ErrorCode!=NO_ERROR) {
72 Finley_Mesh_dealloc(out);
73 return NULL;
74 }
75
76
77 /* allocate tables: */
78
79 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
80 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
81 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
82 if (Finley_ErrorCode!=NO_ERROR) {
83 Finley_Mesh_dealloc(out);
84 return NULL;
85 }
86
87 /* set nodes: */
88
89 #pragma omp parallel for private(i0,i1,i2,k)
90 for (i2=0;i2<N2;i2++) {
91 for (i1=0;i1<N1;i1++) {
92 for (i0=0;i0<N0;i0++) {
93 k=i0+N0*i1+N0*N1*i2;
94 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
95 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
96 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
97 out->Nodes->Id[k]=k;
98 out->Nodes->Tag[k]=0;
99 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1) +N0*N1*(i2%NDOF2);
100 }
101 }
102 }
103 /* tags for the faces: */
104 if (!periodic[2]) {
105 for (i1=0;i1<N1;i1++) {
106 for (i0=0;i0<N0;i0++) {
107 out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;
108 out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;
109 }
110 }
111 }
112 if (!periodic[1]) {
113 for (i2=0;i2<N2;i2++) {
114 for (i0=0;i0<N0;i0++) {
115 out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;
116 out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;
117 }
118 }
119 }
120 if (!periodic[0]) {
121 for (i2=0;i2<N2;i2++) {
122 for (i1=0;i1<N1;i1++) {
123 out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;
124 out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;
125 }
126 }
127 }
128
129 /* set the elements: */
130
131 #pragma omp parallel for private(i0,i1,i2,k,node0)
132 for (i2=0;i2<NE2;i2++) {
133 for (i1=0;i1<NE1;i1++) {
134 for (i0=0;i0<NE0;i0++) {
135 k=i0+NE0*i1+NE0*NE1*i2;
136 node0=i0+i1*N0+N0*N1*i2;
137
138 out->Elements->Id[k]=k;
139 out->Elements->Tag[k]=0;
140 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
141
142 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
143 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
144 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
145 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
146 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
147 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
148 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
149 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
150
151 }
152 }
153 }
154 out->Elements->numColors=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0)+1;
155
156 /* face elements: */
157
158 if (useElementsOnFace) {
159 NUMNODES=8;
160 } else {
161 NUMNODES=4;
162 }
163 totalNECount=NE0*NE1*NE2;
164 faceNECount=0;
165
166 /* these are the quadrilateral elements on boundary 1 (x3=0): */
167 if (!periodic[2]) {
168 /* ** elements on boundary 100 (x3=0): */
169
170 #pragma omp parallel for private(i0,i1,k,node0)
171 for (i1=0;i1<NE1;i1++) {
172 for (i0=0;i0<NE0;i0++) {
173 k=i0+NE0*i1+faceNECount;
174 node0=i0+i1*N0;
175
176 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
177 out->FaceElements->Tag[k]=100;
178 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
179
180 if (useElementsOnFace) {
181 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
182 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
183 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
184 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
185 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
186 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
187 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
188 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
189 } else {
190 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
191 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
192 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
193 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
194 }
195 }
196 }
197 totalNECount+=NE1*NE0;
198 faceNECount+=NE1*NE0;
199
200 /* ** elements on boundary 200 (x3=1) */
201
202 #pragma omp parallel for private(i0,i1,k,node0)
203 for (i1=0;i1<NE1;i1++) {
204 for (i0=0;i0<NE0;i0++) {
205 k=i0+NE0*i1+faceNECount;
206 node0=i0+i1*N0+N0*N1*(NE2-1);
207
208 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
209 out->FaceElements->Tag[k]=200;
210 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
211
212 if (useElementsOnFace) {
213 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
214 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
215 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
216 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
217 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
218 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
219 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
220 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
221 } else {
222 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
223 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
224 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
225 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
226 }
227
228
229 }
230 }
231 totalNECount+=NE1*NE0;
232 faceNECount+=NE1*NE0;
233 }
234 if (!periodic[0]) {
235 /* ** elements on boundary 001 (x1=0): */
236
237 #pragma omp parallel for private(i1,i2,k,node0)
238 for (i2=0;i2<NE2;i2++) {
239 for (i1=0;i1<NE1;i1++) {
240 k=i1+NE1*i2+faceNECount;
241 node0=i1*N0+N0*N1*i2;
242
243 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
244 out->FaceElements->Tag[k]=1;
245 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
246
247 if (useElementsOnFace) {
248 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
249 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
250 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
251 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
252 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
253 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
254 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
255 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
256 } else {
257 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
258 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
259 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
260 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
261 }
262 }
263 }
264 totalNECount+=NE1*NE2;
265 faceNECount+=NE1*NE2;
266
267 /* ** elements on boundary 002 (x1=1): */
268
269 #pragma omp parallel for private(i1,i2,k,node0)
270 for (i2=0;i2<NE2;i2++) {
271 for (i1=0;i1<NE1;i1++) {
272 k=i1+NE1*i2+faceNECount;
273 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
274
275 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
276 out->FaceElements->Tag[k]=2;
277 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
278
279 if (useElementsOnFace) {
280 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
281 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
282 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
283 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
284 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
285 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
286 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
287 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
288 } else {
289 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
290 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
291 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
292 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
293 }
294 }
295 }
296 totalNECount+=NE1*NE2;
297 faceNECount+=NE1*NE2;
298 }
299 if (!periodic[1]) {
300 /* ** elements on boundary 010 (x2=0): */
301
302 #pragma omp parallel for private(i0,i2,k,node0)
303 for (i2=0;i2<NE2;i2++) {
304 for (i0=0;i0<NE0;i0++) {
305 k=i0+NE0*i2+faceNECount;
306 node0=i0+N0*N1*i2;
307
308 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
309 out->FaceElements->Tag[k]=10;
310 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
311
312 if (useElementsOnFace) {
313 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
314 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
315 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
316 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
317 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
318 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
319 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
320 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
321 } else {
322 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
323 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
324 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
325 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
326 }
327 }
328 }
329 totalNECount+=NE0*NE2;
330 faceNECount+=NE0*NE2;
331
332 /* ** elements on boundary 020 (x2=1): */
333
334 #pragma omp parallel for private(i0,i2,k,node0)
335 for (i2=0;i2<NE2;i2++) {
336 for (i0=0;i0<NE0;i0++) {
337 k=i0+NE0*i2+faceNECount;
338 node0=i0+(NE1-1)*N0+N0*N1*i2;
339
340 out->FaceElements->Tag[k]=20;
341 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
342 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
343
344 if (useElementsOnFace) {
345 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
346 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
347 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
348 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
349 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
350 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
351 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
352 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
353 } else {
354 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
355 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
356 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
357 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
358 }
359
360 }
361 }
362 totalNECount+=NE0*NE2;
363 faceNECount+=NE0*NE2;
364 }
365 out->FaceElements->numColors=24;
366
367 /* face elements done: */
368
369 /* condense the nodes: */
370
371 Finley_Mesh_resolveNodeIds(out);
372
373 /* prepare mesh for further calculatuions:*/
374 Finley_Mesh_prepare(out) ;
375
376 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
377
378 if (Finley_ErrorCode!=NO_ERROR) {
379 Finley_Mesh_dealloc(out);
380 return NULL;
381 }
382 return out;
383 }
384
385 /*
386 * $Log$
387 * Revision 1.1 2004/10/26 06:53:57 jgs
388 * Initial revision
389 *
390 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
391 * Initial version of eys using boost-python.
392 *
393 *
394 */
395

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26