/[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 123 - (show annotations)
Fri Jul 8 04:08:13 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/finley/src/finleyC/Mesh_hex8.c
File MIME type: text/plain
File size: 14188 byte(s)
Merge of development branch back to main trunk on 2005-07-08

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26