/[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 964 - (show annotations)
Tue Feb 13 05:10:26 2007 UTC (12 years, 7 months ago) by gross
File MIME type: text/plain
File size: 37555 byte(s)
The set/getRefVal functions of Data objects have been removed (mainly to avoid later problems with MPI).
Moreover, a faster access to the reference id of samples has been introduced. I don't think that anybody will
profit form this at this stage but it will allow a faster dump of data objects.


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 #ifdef PASO_MPI
33 /* get the number of nodes/elements for domain with rank=rank, of size processors
34 where n is the total number of nodes/elements in the global domain */
35 static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36 {
37 rank = size-rank-1;
38
39 if( rank < n%size )
40 return (index_t)floor(n/size)+1;
41 return (index_t)floor(n/size);
42 }
43
44
45 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46 /* A bit messy, but it only has to be done once... */
47 static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal )
48 {
49 index_t i0;
50 dim_t numNodesGlobal = numElementsGlobal+1;
51
52 (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53
54 numElementsLocal[0] = numNodesLocal[0]+1;
55 periodicLocal[0] = periodicLocal[1] = FALSE;
56 nodesExternal[0] = nodesExternal[1] = 1;
57 if( periodic )
58 {
59 if( size==1 )
60 {
61 numElementsLocal[0] = numElementsGlobal;
62 nodesExternal[0] = nodesExternal[1] = 0;
63 periodicLocal[0] = periodicLocal[1] = TRUE;
64 }
65 else
66 {
67 if( rank==0 )
68 {
69 periodicLocal[0] = TRUE;
70 numNodesLocal[0]++;
71 }
72 if( rank==(size-1) )
73 {
74 periodicLocal[1] = TRUE;
75 numNodesLocal[0]--;
76 numElementsLocal[0]--;
77 }
78 }
79 }
80 else if( !periodic )
81 {
82 if( rank==0 ){
83 nodesExternal[0]--;
84 numElementsLocal[0]--;
85 }
86 if( rank==(size-1) )
87 {
88 nodesExternal[1]--;
89 numElementsLocal[0]--;
90 }
91 }
92 numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93
94 numElementsInternal[0] = numElementsLocal[0];
95 if( (rank==0) && (rank==size-1) );
96
97 else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98 numElementsInternal[0] -= 1;
99 else
100 numElementsInternal[0] -= 2;
101
102 firstNode[0] = 0;
103 for( i0=0; i0<rank; i0++ )
104 firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105
106 numDOFLocal[0] = numNodesLocal[0];
107 if( periodicLocal[0] )
108 {
109 numDOFLocal[0]--;
110 }
111 DOFExternal[0] = nodesExternal[0];
112 DOFExternal[1] = nodesExternal[1];
113 }
114
115 #endif
116
117 #ifdef PASO_MPI
118 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119 #else
120 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
121 #endif
122 {
123 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
124 index_t node0;
125 Finley_Mesh* out;
126 char name[50];
127 double time0=Finley_timer();
128 NE0=MAX(1,numElements[0]);
129 NE1=MAX(1,numElements[1]);
130 NE2=MAX(1,numElements[2]);
131 N0=NE0+1;
132 N1=NE1+1;
133 N2=NE2+1;
134
135 if (N0<=MIN(N1,N2)) {
136 if (N1 <= N2) {
137 M0=1;
138 M1=N0;
139 M2=N0*N1;
140 } else {
141 M0=1;
142 M2=N0;
143 M1=N0*N2;
144 }
145 } else if (N1<=MIN(N2,N0)) {
146 if (N2 <= N0) {
147 M1=1;
148 M2=N1;
149 M0=N2*N1;
150 } else {
151 M1=1;
152 M0=N1;
153 M2=N1*N0;
154 }
155 } else {
156 if (N0 <= N1) {
157 M2=1;
158 M0=N2;
159 M1=N2*N0;
160 } else {
161 M2=1;
162 M1=N2;
163 M0=N1*N2;
164 }
165 }
166
167
168 NFaceElements=0;
169 if (!periodic[0]) {
170 NDOF0=N0;
171 NFaceElements+=2*NE1*NE2;
172 } else {
173 NDOF0=N0-1;
174 }
175 if (!periodic[1]) {
176 NDOF1=N1;
177 NFaceElements+=2*NE0*NE2;
178 } else {
179 NDOF1=N1-1;
180 }
181 if (!periodic[2]) {
182 NDOF2=N2;
183 NFaceElements+=2*NE0*NE1;
184 } else {
185 NDOF2=N2-1;
186 }
187
188 /* allocate mesh: */
189
190 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
191
192 #ifndef PASO_MPI
193 out=Finley_Mesh_alloc(name,3,order);
194 #else
195 out=Finley_Mesh_alloc(name,3,order,mpi_info);
196 #endif
197 if (! Finley_noError()) return NULL;
198
199 #ifdef PASO_MPI
200 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
201 if (useElementsOnFace) {
202 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
203 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
204 } else {
205 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
206 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
207 }
208 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
209 #else
210 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
211 if (useElementsOnFace) {
212 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
213 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
214 } else {
215 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
216 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
217 }
218 out->Points=Finley_ElementFile_alloc(Point1,out->order);
219 #endif
220 if (! Finley_noError()) {
221 Finley_Mesh_dealloc(out);
222 return NULL;
223 }
224
225
226 /* allocate tables: */
227 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
228 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
229 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
230 #ifdef PASO_MPI
231 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
232 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1*NE2, NE0*NE1*NE2);
233 Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );
234 #endif
235 if (! Finley_noError()) {
236 Finley_Mesh_dealloc(out);
237 return NULL;
238 }
239
240 /* set nodes: */
241
242 #pragma omp parallel for private(i0,i1,i2,k)
243 for (i2=0;i2<N2;i2++) {
244 for (i1=0;i1<N1;i1++) {
245 for (i0=0;i0<N0;i0++) {
246 k=M0*i0+M1*i1+M2*i2;
247 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
248 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
249 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
250 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
251 out->Nodes->Tag[k]=0;
252 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
253 #ifdef PASO_MPI
254 out->Nodes->Dom[k]=NODE_INTERNAL;
255 #endif
256 }
257 }
258 }
259 /* tags for the faces: */
260 /* tags for the faces: */
261 if (!periodic[2]) {
262 for (i1=0;i1<N1;i1++) {
263 for (i0=0;i0<N0;i0++) {
264 out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
265 out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
266 }
267 }
268 }
269 if (!periodic[1]) {
270 for (i2=0;i2<N2;i2++) {
271 for (i0=0;i0<N0;i0++) {
272 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
273 out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
274 }
275 }
276 }
277 if (!periodic[0]) {
278 for (i2=0;i2<N2;i2++) {
279 for (i1=0;i1<N1;i1++) {
280 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
281 out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
282 }
283 }
284 }
285
286 /* set the elements: */
287
288 #pragma omp parallel for private(i0,i1,i2,k,node0)
289 for (i2=0;i2<NE2;i2++) {
290 for (i1=0;i1<NE1;i1++) {
291 for (i0=0;i0<NE0;i0++) {
292 k=i0+NE0*i1+NE0*NE1*i2;
293 node0=i0+i1*N0+N0*N1*i2;
294
295 out->Elements->Id[k]=k;
296 out->Elements->Tag[k]=0;
297 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
298 #ifdef PASO_MPI
299 out->Elements->Dom[k]=ELEMENT_INTERNAL;
300 #endif
301
302 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
303 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
304 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
305 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
306 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
307 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
308 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
309 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
310
311 }
312 }
313 }
314 out->Elements->minColor=0;
315 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
316
317 /* face elements: */
318
319 if (useElementsOnFace) {
320 NUMNODES=8;
321 } else {
322 NUMNODES=4;
323 }
324 totalNECount=NE0*NE1*NE2;
325 faceNECount=0;
326
327 /* these are the quadrilateral elements on boundary 1 (x3=0): */
328 if (!periodic[2]) {
329 /* ** elements on boundary 100 (x3=0): */
330
331 #pragma omp parallel for private(i0,i1,k,node0)
332 for (i1=0;i1<NE1;i1++) {
333 for (i0=0;i0<NE0;i0++) {
334 k=i0+NE0*i1+faceNECount;
335 node0=i0+i1*N0;
336
337 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
338 out->FaceElements->Tag[k]=100;
339 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
340 #ifdef PASO_MPI
341 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
342 #endif
343
344 if (useElementsOnFace) {
345 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
346 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
347 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
348 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
349 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
350 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
351 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
352 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
353 } else {
354 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
355 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
356 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
357 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
358 }
359 }
360 }
361 totalNECount+=NE1*NE0;
362 faceNECount+=NE1*NE0;
363
364 /* ** elements on boundary 200 (x3=1) */
365
366 #pragma omp parallel for private(i0,i1,k,node0)
367 for (i1=0;i1<NE1;i1++) {
368 for (i0=0;i0<NE0;i0++) {
369 k=i0+NE0*i1+faceNECount;
370 node0=i0+i1*N0+N0*N1*(NE2-1);
371
372 out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
373 out->FaceElements->Tag[k]=200;
374 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
375 #ifdef PASO_MPI
376 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
377 #endif
378
379 if (useElementsOnFace) {
380 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
381 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
382 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
383 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
384 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
385 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
386 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
387 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
388 } else {
389 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
390 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
391 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
392 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
393 }
394
395
396 }
397 }
398 totalNECount+=NE1*NE0;
399 faceNECount+=NE1*NE0;
400 }
401 if (!periodic[0]) {
402 /* ** elements on boundary 001 (x1=0): */
403
404 #pragma omp parallel for private(i1,i2,k,node0)
405 for (i2=0;i2<NE2;i2++) {
406 for (i1=0;i1<NE1;i1++) {
407 k=i1+NE1*i2+faceNECount;
408 node0=i1*N0+N0*N1*i2;
409
410 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
411 out->FaceElements->Tag[k]=1;
412 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
413 #ifdef PASO_MPI
414 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
415 #endif
416
417 if (useElementsOnFace) {
418 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
419 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
420 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
421 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
422 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
423 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
424 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
425 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
426 } else {
427 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
428 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
429 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
430 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
431 }
432 }
433 }
434 totalNECount+=NE1*NE2;
435 faceNECount+=NE1*NE2;
436
437 /* ** elements on boundary 002 (x1=1): */
438
439 #pragma omp parallel for private(i1,i2,k,node0)
440 for (i2=0;i2<NE2;i2++) {
441 for (i1=0;i1<NE1;i1++) {
442 k=i1+NE1*i2+faceNECount;
443 node0=(NE0-1)+i1*N0+N0*N1*i2 ;
444
445 out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
446 out->FaceElements->Tag[k]=2;
447 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
448 #ifdef PASO_MPI
449 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
450 #endif
451
452 if (useElementsOnFace) {
453 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
454 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
455 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
456 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
457 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
458 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
459 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
460 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
461 } else {
462 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
463 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
464 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
465 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
466 }
467 }
468 }
469 totalNECount+=NE1*NE2;
470 faceNECount+=NE1*NE2;
471 }
472 if (!periodic[1]) {
473 /* ** elements on boundary 010 (x2=0): */
474
475 #pragma omp parallel for private(i0,i2,k,node0)
476 for (i2=0;i2<NE2;i2++) {
477 for (i0=0;i0<NE0;i0++) {
478 k=i0+NE0*i2+faceNECount;
479 node0=i0+N0*N1*i2;
480
481 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
482 out->FaceElements->Tag[k]=10;
483 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
484 #ifdef PASO_MPI
485 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
486 #endif
487
488 if (useElementsOnFace) {
489 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
490 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
491 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
492 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
493 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
494 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
495 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
496 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
497 } else {
498 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
499 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
500 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
501 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
502 }
503 }
504 }
505 totalNECount+=NE0*NE2;
506 faceNECount+=NE0*NE2;
507
508 /* ** elements on boundary 020 (x2=1): */
509
510 #pragma omp parallel for private(i0,i2,k,node0)
511 for (i2=0;i2<NE2;i2++) {
512 for (i0=0;i0<NE0;i0++) {
513 k=i0+NE0*i2+faceNECount;
514 node0=i0+(NE1-1)*N0+N0*N1*i2;
515
516 out->FaceElements->Tag[k]=20;
517 out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
518 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
519 #ifdef PASO_MPI
520 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
521 #endif
522
523 if (useElementsOnFace) {
524 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
525 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
526 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
527 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
528 out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
529 out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
530 out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
531 out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
532 } else {
533 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
534 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
535 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
536 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
537 }
538
539 }
540 }
541 totalNECount+=NE0*NE2;
542 faceNECount+=NE0*NE2;
543 }
544 out->FaceElements->minColor=0;
545 out->FaceElements->maxColor=23;
546
547 #ifdef PASO_MPI
548 Finley_ElementFile_setDomainFlags( out->Elements );
549 Finley_ElementFile_setDomainFlags( out->FaceElements );
550 Finley_ElementFile_setDomainFlags( out->ContactElements );
551 Finley_ElementFile_setDomainFlags( out->Points );
552
553 /* reorder the degrees of freedom */
554 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
555 #endif
556
557 /* condense the nodes: */
558 Finley_Mesh_resolveNodeIds(out);
559
560 /* prepare mesh for further calculatuions:*/
561 Finley_Mesh_prepare(out) ;
562
563 #ifdef Finley_TRACE
564 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
565 #endif
566 if (Finley_noError()) {
567 if ( ! Finley_Mesh_isPrepared(out)) {
568 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
569 }
570 } else {
571 Finley_Mesh_dealloc(out);
572 }
573 return out;
574 }
575
576 #ifdef PASO_MPI
577 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
578 {
579 dim_t N0,N1,N2,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
580 dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
581 bool_t dom_left, dom_right, dom_internal;
582 index_t firstNode=0, DOFcount=0, node0, node1, node2;
583 index_t targetDomain=-1, firstNodeConstruct, j;
584 bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
585 index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
586 Finley_Mesh* out;
587
588 char name[50];
589 Paso_MPIInfo *mpi_info = NULL;
590 double time0=Finley_timer();
591
592 index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
593 index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
594 index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
595 index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
596 index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
597 index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
598 NE0=MAX(1,numElements[0]);
599 NE1=MAX(1,numElements[1]);
600 NE2=MAX(1,numElements[2]);
601 N0=NE0+1;
602 N1=NE1+1;
603 N2=NE2+1;
604
605
606 /* get MPI information */
607 mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
608 if (! Finley_noError())
609 return NULL;
610
611 /* use the serial version to generate the mesh for the 1-CPU case */
612 if( mpi_info->size==1 )
613 {
614 out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info );
615 return out;
616 }
617
618 if( mpi_info->rank==0 )
619 domLeft = TRUE;
620 if( mpi_info->rank==mpi_info->size-1 )
621 domRight = TRUE;
622 if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
623 domInternal = TRUE;
624
625 /* dimensions of the local subdomain */
626 domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
627
628 /* count Degrees of Freedom along each axis */
629 NDOF0 = N0 - periodic[0];
630 NDOF1 = N1 - periodic[1];
631 NDOF2 = N2 - periodic[2];
632
633 /* count face elements */
634 /* internal face elements */
635 NFaceElements = 0;
636 if( !periodic[0] )
637 NFaceElements += (domLeft+domRight)*NE1*NE2;
638 if( !periodic[1] )
639 NFaceElements += 2*numElementsLocal*NE2;
640 if( !periodic[2] )
641 NFaceElements += 2*numElementsLocal*NE1;
642
643 boundaryLeft = !domLeft || periodicLocal[0];
644 boundaryRight = !domRight || periodicLocal[1];
645 N0t = numNodesLocal + boundaryRight + boundaryLeft;
646 NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
647 firstNodeConstruct = firstNode - boundaryLeft;
648 firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
649
650 /* allocate mesh: */
651 sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
652
653 out=Finley_Mesh_alloc(name,3,order,mpi_info);
654 if (! Finley_noError()) return NULL;
655
656 out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
657 if (useElementsOnFace) {
658 out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
659 out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
660 } else {
661 out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
662 out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
663 }
664 out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
665
666 if (! Finley_noError()) {
667 Finley_Mesh_dealloc(out);
668 return NULL;
669 }
670
671 /* allocate tables: */
672 Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
673 Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
674 Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
675
676 Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
677 Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );
678 Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );
679 if (! Finley_noError()) {
680 Finley_Mesh_dealloc(out);
681 return NULL;
682 }
683
684 /* set nodes: */
685 /* INTERNAL/BOUNDARY NODES */
686 k=0;
687 #pragma omp parallel for private(i0,i1,i2,k)
688 for (i2=0;i2<N2;i2++) {
689 for (i1=0;i1<N1;i1++) {
690 for (i0=0;i0<N0t;i0++,k++) {
691 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
692 out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
693 out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
694 out->Nodes->Id[k]=k;
695 out->Nodes->Tag[k]=0;
696 out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
697 out->Nodes->Dom[k]=NODE_INTERNAL;
698 }
699 }
700 }
701
702 /* mark the nodes that reference external and boundary DOF as such */
703 if( boundaryLeft ){
704 for( i1=0; i1<N1; i1++ )
705 for( i2=0; i2<N2; i2++ ) {
706 out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
707 out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
708 }
709 }
710 if( boundaryRight ){
711 for( i1=0; i1<N1; i1++ )
712 for( i2=0; i2<N2; i2++ ) {
713 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
714 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
715 }
716 }
717 if( periodicLocal[0] ){
718 for( i1=0; i1<N1; i1++ )
719 for( i2=0; i2<N2; i2++ ) {
720 out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
721 out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
722 }
723 }
724
725 /* tag Nodes that are referenced by face elements */
726 if (!periodic[2]) {
727 for (i1=0;i1<N1;i1++) {
728 for (i0=0;i0<N0t;i0++) {
729 out->Nodes->Tag[i0 + N0t*i1]+=100;
730 out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
731 }
732 }
733 }
734 if (!periodic[1]) {
735 for (i2=0;i2<N2;i2++) {
736 for (i0=0;i0<N0t;i0++) {
737 out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
738 out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
739 }
740 }
741 }
742 if (!periodic[0] && !domInternal ) {
743 for (i2=0;i2<N2;i2++) {
744 for (i1=0;i1<N1;i1++) {
745 if( domLeft )
746 out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
747 if( domRight )
748 out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
749 }
750 }
751 }
752
753 /* form the boudary communication information */
754 forwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
755 backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
756 if( !(mpi_info->size==2 && periodicLocal[0])){
757 if( boundaryLeft ) {
758 targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
759 for( i2=0; i2<NDOF2; i2++ ){
760 for( i1=0; i1<NDOF1; i1++ ){
761 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
762 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
763 }
764 }
765 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
766 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
767 }
768 if( boundaryRight ) {
769 targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
770 for( i2=0; i2<NDOF2; i2++ ){
771 for( i1=0; i1<NDOF1; i1++ ){
772 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
773 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
774 }
775 }
776 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
777 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
778 }
779 } else{
780 /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
781 targetDomain = 1;
782
783 for( i2=0; i2<NDOF2; i2++ ){
784 for( i1=0; i1<NDOF1; i1++ ){
785 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
786 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
787 }
788 }
789 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
790 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
791
792 for( i2=0; i2<NDOF2; i2++ ){
793 for( i1=0; i1<NDOF1; i1++ ){
794 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
795 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
796 }
797 }
798 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
799 Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
800 }
801 MEMFREE( forwardDOF );
802 MEMFREE( backwardDOF );
803 /* set the elements: */
804
805 /* INTERNAL elements */
806 k = 0;
807 #pragma omp parallel for private(i0,i1,i2,k,node0)
808 for (i2=0;i2<NE2;i2++) {
809 for (i1=0;i1<NE1;i1++) {
810 for (i0=0;i0<numElementsLocal;i0++,k++) {
811 node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t : i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
812
813 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
814 out->Elements->Tag[k]=0;
815 out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
816 out->Elements->Dom[k]=ELEMENT_INTERNAL;
817
818 out->Elements->Nodes[INDEX2(0,k,8)]=node0;
819 out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
820 out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
821 out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
822 out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
823 out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
824 out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
825 out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
826
827 }
828 }
829 }
830 out->Elements->minColor=0;
831 out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
832 if( boundaryLeft )
833 for( i2=0; i2<NE2; i2++ )
834 for( i1=0; i1<NE1; i1++ )
835 out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
836 if( boundaryRight )
837 for( i2=0; i2<NE2; i2++ )
838 for( i1=0; i1<NE1; i1++ )
839 out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
840
841 Finley_ElementFile_setDomainFlags( out->Elements );
842
843 /* face elements: */
844 if (useElementsOnFace) {
845 NUMNODES=8;
846 } else {
847 NUMNODES=4;
848 }
849 totalNECount=out->Elements->numElements;
850 faceNECount=0;
851 idCount = totalNECount;
852
853 /* these are the quadrilateral elements on boundary 1 (x3=0): */
854 numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
855 if (!periodic[2]) {
856 /* elements on boundary 100 (x3=0): */
857
858 #pragma omp parallel for private(i0,i1,k)
859 for (i1=0;i1<NE1;i1++) {
860 for (i0=0; i0<numElementsLocal; i0++) {
861 k=i0+numElementsLocal*i1+faceNECount;
862 kk=i0 + i1*numElementsLocal;
863 facePerm = face2;
864
865 out->FaceElements->Id[k]=idCount++;
866 out->FaceElements->Tag[k]=100;
867 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
868 out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
869
870 for( j=0; j<NUMNODES; j++ )
871 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
872 }
873 }
874 if( boundaryLeft ){
875 for( i1=0; i1<NE1; i1++ )
876 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
877 if( periodicLocal[0] )
878 for( i1=0; i1<NE1; i1++ )
879 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
880 }
881 if( boundaryRight )
882 for( i1=0; i1<NE1; i1++ )
883 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
884 totalNECount+=NE1*numElementsLocal;
885 faceNECount+=NE1*numElementsLocal;
886
887 /* ** elements on boundary 200 (x3=1) */
888
889 #pragma omp parallel for private(i0,i1,k)
890 for (i1=0;i1<NE1;i1++) {
891 for (i0=0;i0<numElementsLocal;i0++) {
892 k=i0+numElementsLocal*i1+faceNECount;
893 kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
894 facePerm = face3;
895
896 out->FaceElements->Id[k]=idCount++;
897 out->FaceElements->Tag[k]=200;
898 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
899 out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
900
901 for( j=0; j<NUMNODES; j++ )
902 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
903 }
904 }
905 if( boundaryLeft ){
906 for( i1=0; i1<NE1; i1++ )
907 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
908 if( periodicLocal[0] )
909 for( i1=0; i1<NE1; i1++ )
910 out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
911 }
912 if( boundaryRight )
913 for( i1=0; i1<NE1; i1++ )
914 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
915 totalNECount+=NE1*numElementsLocal;
916 faceNECount+=NE1*numElementsLocal;
917 }
918 if (!periodic[0] && !domInternal) {
919 /* ** elements on boundary 001 (x1=0): */
920 if( domLeft ){
921 #pragma omp parallel for private(i1,i2,k)
922 for (i2=0;i2<NE2;i2++) {
923 for (i1=0;i1<NE1;i1++) {
924 k=i1+NE1*i2+faceNECount;
925 kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
926 facePerm = face0;
927
928 out->FaceElements->Id[k]=idCount++;
929 out->FaceElements->Tag[k]=1;
930 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
931 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
932
933 for( j=0; j<NUMNODES; j++ )
934 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
935 }
936 }
937 totalNECount+=NE1*NE2;
938 faceNECount+=NE1*NE2;
939 }
940 /* ** elements on boundary 002 (x1=1): */
941 if( domRight ) {
942 #pragma omp parallel for private(i1,i2,k)
943 for (i2=0;i2<NE2;i2++) {
944 for (i1=0;i1<NE1;i1++) {
945 k=i1+NE1*i2+faceNECount;
946 kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
947 facePerm = face1;
948
949 out->FaceElements->Id[k]=idCount++;
950 out->FaceElements->Tag[k]=2;
951 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
952 out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
953
954 for( j=0; j<NUMNODES; j++ )
955 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
956 }
957 }
958 totalNECount+=NE1*NE2;
959 faceNECount+=NE1*NE2;
960 }
961 }
962 if (!periodic[1]) {
963 /* ** elements on boundary 010 (x2=0): */
964
965 #pragma omp parallel for private(i0,i2,k)
966 for (i2=0;i2<NE2;i2++) {
967 for (i0=0;i0<numElementsLocal;i0++) {
968 k=i0+numElementsLocal*i2+faceNECount;
969 kk=i0+numElementsLocal*NE1*i2;
970 facePerm = face4;
971
972 out->FaceElements->Id[k]=idCount++;
973 out->FaceElements->Tag[k]=10;
974 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
975 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
976
977 for( j=0; j<NUMNODES; j++ )
978 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
979 }
980 }
981 if( boundaryLeft ){
982 for( i2=0; i2<NE2; i2++ )
983 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
984 if( periodicLocal[0] )
985 for( i2=0; i2<NE2; i2++ )
986 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
987 }
988 if( boundaryRight )
989 for( i2=0; i2<NE2; i2++ )
990 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
991 totalNECount+=numElementsLocal*NE2;
992 faceNECount+=numElementsLocal*NE2;
993
994 /* ** elements on boundary 020 (x2=1): */
995
996 #pragma omp parallel for private(i0,i2,k)
997 for (i2=0;i2<NE2;i2++) {
998 for (i0=0;i0<numElementsLocal;i0++) {
999 k=i0+numElementsLocal*i2+faceNECount;
1000 kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
1001 facePerm = face5;
1002
1003 out->FaceElements->Tag[k]=20;
1004 out->FaceElements->Id[k]=idCount++;
1005 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1006 out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1007
1008 for( j=0; j<NUMNODES; j++ )
1009 out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
1010 }
1011 }
1012 if( boundaryLeft ){
1013 for( i2=0; i2<NE2; i2++ )
1014 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1015 if( periodicLocal[0] )
1016 for( i2=0; i2<NE2; i2++ )
1017 out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1018 }
1019 if( boundaryRight )
1020 for( i2=0; i2<NE2; i2++ )
1021 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1022 totalNECount+=numElementsLocal*NE2;
1023 faceNECount+=numElementsLocal*NE2;
1024 }
1025 out->FaceElements->elementDistribution->numInternal = faceNECount;
1026
1027 out->FaceElements->minColor=0;
1028 out->FaceElements->maxColor=23;
1029 out->FaceElements->numElements=faceNECount;
1030
1031 Finley_ElementFile_setDomainFlags( out->FaceElements );
1032
1033 /* setup distribution info for other elements */
1034 Finley_ElementFile_setDomainFlags( out->ContactElements );
1035 Finley_ElementFile_setDomainFlags( out->Points );
1036
1037 /* reorder the degrees of freedom */
1038 Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1039
1040 /* condense the nodes: */
1041 Finley_Mesh_resolveNodeIds(out);
1042 if( !Finley_MPI_noError(mpi_info) )
1043 {
1044 Paso_MPIInfo_dealloc( mpi_info );
1045 Finley_Mesh_dealloc(out);
1046 return NULL;
1047 }
1048
1049 /* prepare mesh for further calculatuions:*/
1050 Finley_Mesh_prepare(out);
1051 if( !Finley_MPI_noError(mpi_info) )
1052 {
1053 Paso_MPIInfo_dealloc( mpi_info );
1054 Finley_Mesh_dealloc(out);
1055 return NULL;
1056 }
1057
1058 /* free up memory */
1059 Paso_MPIInfo_dealloc( mpi_info );
1060
1061 #ifdef Finley_TRACE
1062 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1063 #endif
1064 if (Finley_noError()) {
1065 if ( ! Finley_Mesh_isPrepared(out) ) {
1066 Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
1067 }
1068 }
1069 return out;
1070 }
1071 #endif
1072

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26