/[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 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (13 years, 11 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex8.c
File MIME type: text/plain
File size: 15338 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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;
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 NFaceElements=0;
48 if (!periodic[0]) {
49 NDOF0=N0;
50 NFaceElements+=2*NE1*NE2;
51 } else {
52 NDOF0=N0-1;
53 }
54 if (!periodic[1]) {
55 NDOF1=N1;
56 NFaceElements+=2*NE0*NE2;
57 } else {
58 NDOF1=N1-1;
59 }
60 if (!periodic[2]) {
61 NDOF2=N2;
62 NFaceElements+=2*NE0*NE1;
63 } else {
64 NDOF2=N2-1;
65 }
66
67 /* allocate mesh: */
68
69 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
70 out=Finley_Mesh_alloc(name,3,order);
71 if (! Finley_noError()) return NULL;
72
73 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
74 if (useElementsOnFace) {
75 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
76 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
77 } else {
78 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
79 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
80 }
81 out->Points=Finley_ElementFile_alloc(Point1,out->order);
82 if (! Finley_noError()) {
83 Finley_Mesh_dealloc(out);
84 return NULL;
85 }
86
87
88 /* allocate tables: */
89
90 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
91 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
92 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
93 if (! Finley_noError()) {
94 Finley_Mesh_dealloc(out);
95 return NULL;
96 }
97
98 /* set nodes: */
99
100 #pragma omp parallel for private(i0,i1,i2,k)
101 for (i2=0;i2<N2;i2++) {
102 for (i1=0;i1<N1;i1++) {
103 for (i0=0;i0<N0;i0++) {
104 k=i0+N0*i1+N0*N1*i2;
105 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
106 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
107 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
108 out->Nodes->Id[k]=k;
109 out->Nodes->Tag[k]=0;
110 out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1) +N0*N1*(i2%NDOF2);
111 }
112 }
113 }
114 /* tags for the faces: */
115 if (!periodic[2]) {
116 for (i1=0;i1<N1;i1++) {
117 for (i0=0;i0<N0;i0++) {
118 out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;
119 out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;
120 }
121 }
122 }
123 if (!periodic[1]) {
124 for (i2=0;i2<N2;i2++) {
125 for (i0=0;i0<N0;i0++) {
126 out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;
127 out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;
128 }
129 }
130 }
131 if (!periodic[0]) {
132 for (i2=0;i2<N2;i2++) {
133 for (i1=0;i1<N1;i1++) {
134 out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;
135 out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;
136 }
137 }
138 }
139
140 /* set the elements: */
141
142 #pragma omp parallel for private(i0,i1,i2,k,node0)
143 for (i2=0;i2<NE2;i2++) {
144 for (i1=0;i1<NE1;i1++) {
145 for (i0=0;i0<NE0;i0++) {
146 k=i0+NE0*i1+NE0*NE1*i2;
147 node0=i0+i1*N0+N0*N1*i2;
148
149 out->Elements->Id[k]=k;
150 out->Elements->Tag[k]=0;
151 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
152
153 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
154 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
155 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
156 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
157 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
158 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
159 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
160 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
161
162 }
163 }
164 }
165 out->Elements->minColor=0;
166 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
167
168 /* face elements: */
169
170 if (useElementsOnFace) {
171 NUMNODES=8;
172 } else {
173 NUMNODES=4;
174 }
175 totalNECount=NE0*NE1*NE2;
176 faceNECount=0;
177
178 /* these are the quadrilateral elements on boundary 1 (x3=0): */
179 if (!periodic[2]) {
180 /* ** elements on boundary 100 (x3=0): */
181
182 #pragma omp parallel for private(i0,i1,k,node0)
183 for (i1=0;i1<NE1;i1++) {
184 for (i0=0;i0<NE0;i0++) {
185 k=i0+NE0*i1+faceNECount;
186 node0=i0+i1*N0;
187
188 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
189 out->FaceElements->Tag[k]=100;
190 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
191
192 if (useElementsOnFace) {
193 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
194 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
195 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
196 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
197 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
198 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
199 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
200 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
201 } else {
202 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
203 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
204 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
205 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
206 }
207 }
208 }
209 totalNECount+=NE1*NE0;
210 faceNECount+=NE1*NE0;
211
212 /* ** elements on boundary 200 (x3=1) */
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+N0*N1*(NE2-1);
219
220 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
221 out->FaceElements->Tag[k]=200;
222 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
223
224 if (useElementsOnFace) {
225 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
226 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
227 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
228 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
229 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
230 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
231 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
232 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
233 } else {
234 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
235 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
236 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
237 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
238 }
239
240
241 }
242 }
243 totalNECount+=NE1*NE0;
244 faceNECount+=NE1*NE0;
245 }
246 if (!periodic[0]) {
247 /* ** elements on boundary 001 (x1=0): */
248
249 #pragma omp parallel for private(i1,i2,k,node0)
250 for (i2=0;i2<NE2;i2++) {
251 for (i1=0;i1<NE1;i1++) {
252 k=i1+NE1*i2+faceNECount;
253 node0=i1*N0+N0*N1*i2;
254
255 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
256 out->FaceElements->Tag[k]=1;
257 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
258
259 if (useElementsOnFace) {
260 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
261 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
262 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
263 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
264 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
265 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
266 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
267 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
268 } else {
269 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
270 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
271 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
272 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
273 }
274 }
275 }
276 totalNECount+=NE1*NE2;
277 faceNECount+=NE1*NE2;
278
279 /* ** elements on boundary 002 (x1=1): */
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=(NE0-1)+i1*N0+N0*N1*i2 ;
286
287 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
288 out->FaceElements->Tag[k]=2;
289 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
290
291 if (useElementsOnFace) {
292 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
293 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
294 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
295 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
296 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
297 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
298 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
299 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
300 } else {
301 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
302 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
303 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
304 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
305 }
306 }
307 }
308 totalNECount+=NE1*NE2;
309 faceNECount+=NE1*NE2;
310 }
311 if (!periodic[1]) {
312 /* ** elements on boundary 010 (x2=0): */
313
314 #pragma omp parallel for private(i0,i2,k,node0)
315 for (i2=0;i2<NE2;i2++) {
316 for (i0=0;i0<NE0;i0++) {
317 k=i0+NE0*i2+faceNECount;
318 node0=i0+N0*N1*i2;
319
320 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
321 out->FaceElements->Tag[k]=10;
322 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
323
324 if (useElementsOnFace) {
325 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
326 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
327 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
328 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
329 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
330 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
331 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
332 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
333 } else {
334 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
335 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
336 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
337 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
338 }
339 }
340 }
341 totalNECount+=NE0*NE2;
342 faceNECount+=NE0*NE2;
343
344 /* ** elements on boundary 020 (x2=1): */
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+(NE1-1)*N0+N0*N1*i2;
351
352 out->FaceElements->Tag[k]=20;
353 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
354 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
355
356 if (useElementsOnFace) {
357 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
358 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
359 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
360 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
361 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
362 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
363 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
364 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
365 } else {
366 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
367 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
368 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
369 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
370 }
371
372 }
373 }
374 totalNECount+=NE0*NE2;
375 faceNECount+=NE0*NE2;
376 }
377 out->FaceElements->minColor=0;
378 out->FaceElements->maxColor=23;
379
380 /* face elements done: */
381
382 /* condense the nodes: */
383
384 Finley_Mesh_resolveNodeIds(out);
385
386 /* prepare mesh for further calculatuions:*/
387 Finley_Mesh_prepare(out) ;
388
389 #ifdef Finley_TRACE
390 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
391 #endif
392
393 if (! Finley_noError()) {
394 Finley_Mesh_dealloc(out);
395 return NULL;
396 }
397 return out;
398 }
399
400 /*
401 * Revision 1.3 2005/09/01 03:31:35 jgs
402 * Merge of development branch dev-02 back to main trunk on 2005-09-01
403 *
404 * Revision 1.2.2.2 2005/09/07 06:26:19 gross
405 * the solver from finley are put into the standalone package paso now
406 *
407 * Revision 1.2.2.1 2005/08/24 02:02:18 gross
408 * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.
409 *
410 * Revision 1.2 2005/07/08 04:07:52 jgs
411 * Merge of development branch back to main trunk on 2005-07-08
412 *
413 * Revision 1.1.1.1.2.1 2005/06/29 02:34:52 gross
414 * some changes towards 64 integers in finley
415 *
416 * Revision 1.1.1.1 2004/10/26 06:53:57 jgs
417 * initial import of project esys2
418 *
419 * Revision 1.1.1.1 2004/06/24 04:00:40 johng
420 * Initial version of eys using boost-python.
421 *
422 *
423 */
424

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26