/[escript]/trunk/finley/src/Mesh_hex8.c
ViewVC logotype

Annotation of /trunk/finley/src/Mesh_hex8.c

Parent Directory Parent Directory | Revision Log Revision Log


Revision 782 - (hide annotations)
Tue Jul 18 00:47:47 2006 UTC (13 years, 1 month ago) by bcumming
File MIME type: text/plain
File size: 36621 byte(s)
Large number of changes to Finley for meshing in MPI.

- optimisation and neatening up of rectcanglular mesh generation code
- first and second order 1D, 2D and 3D rectangular meshes are now
  available in finley and escript using MPI.
- reduced meshes now generated in MPI, and interpolation to and from 
  reduced data types now supported.  

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26