/[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 759 - (hide annotations)
Thu Jun 29 01:53:23 2006 UTC (13 years, 3 months ago) by bcumming
File MIME type: text/plain
File size: 64656 byte(s)
- added directory pythonMPI to the source tree. this directory contains
  the c++ wrapper that is used to run python scripts in parallel for the
  MPI version of escript/finley
- updated the SConstruct and ./scons/ess_options.py for conditional MPI
  compilation. To compile the MPI version on ESS uncomment the #define
  PASO_MPI in ./paso/src/Paso.h and add the command line option
  useMPI=yes when running scons.
- fixed a compile time error in the MPI build in  
  finley/src/CPPAdapter/MeshAdapterFactory.cpp 
     

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     void print_mesh_statistics( Finley_Mesh *out )
116     {
117     index_t i, j, N;
118    
119     printf( "\nNodes\n=====\n\n" );
120     printf( "\t%d internal DOF\n\t%d boundary DOF\n\t%d local DOF\n\t%d external DOF\n", out->Nodes->degreeOfFreedomDistribution->numInternal, out->Nodes->degreeOfFreedomDistribution->numBoundary, out->Nodes->degreeOfFreedomDistribution->numLocal, out->Nodes->degreeOfFreedomDistribution->numExternal);
121     for( i=0; i<out->Nodes->numNodes; i++ )
122     printf( "node %d\t: id %d \tDOF %d \t: tag %d \t: coordinates [%3g, %3g, %3g]\n", i, out->Nodes->Id[i], out->Nodes->degreeOfFreedom[i], out->Nodes->Tag[i], out->Nodes->Coordinates[INDEX2(0,i,3)], out->Nodes->Coordinates[INDEX2(1,i,3)], out->Nodes->Coordinates[INDEX2(2,i,3)] );
123    
124     printf( "Elements\n========\n\n" );
125     printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->Elements->elementDistribution->numInternal, out->Elements->elementDistribution->numBoundary, out->Elements->elementDistribution->numLocal );
126     N = out->Elements->ReferenceElement->Type->numNodes;
127     for( i=0; i<out->Elements->numElements; i++ ){
128     printf( "element %d \t: id %d \t: nodes [ %3d, %3d, %3d, %3d, %3d, %3d, %3d, %3d ]", i, out->Elements->Id[i], out->Elements->Nodes[INDEX2(0,i,8)], out->Elements->Nodes[INDEX2(1,i,8)], out->Elements->Nodes[INDEX2(2,i,8)], out->Elements->Nodes[INDEX2(3,i,8)], out->Elements->Nodes[INDEX2(4,i,8)], out->Elements->Nodes[INDEX2(5,i,8)], out->Elements->Nodes[INDEX2(6,i,8)], out->Elements->Nodes[INDEX2(7,i,8)] );
129     printf( " DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(0,i,N)]]] );
130     for( j=1; j<N; j++ )
131     printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(j,i,N)]]] );
132     printf( " ]\n" );
133     }
134    
135     printf( "\nFace Elements\n==============\n\n" );
136     printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->FaceElements->elementDistribution->numInternal, out->FaceElements->elementDistribution->numBoundary, out->FaceElements->elementDistribution->numLocal );
137     N = out->FaceElements->ReferenceElement->Type->numNodes;
138     for( i=0; i<out->FaceElements->numElements; i++ ){
139     printf( "face element %d \t: id %d \t: nodes [ %3d", i, out->FaceElements->Id[i], out->FaceElements->Nodes[INDEX2(0,i,N)] );
140     for( j=1; j<N; j++ )
141     printf( ", %3d", out->FaceElements->Nodes[INDEX2(j,i,N)] );
142     printf( " ] DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(0,i,N)]]] );
143     for( j=1; j<N; j++ )
144     printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,N)]]] );
145     printf( " ]\n" );
146     }
147     }
148    
149     #endif
150    
151     #ifndef PASO_MPI
152     Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
153     #else
154     Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
155     #endif
156     {
157 jgs 153 dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
158 jgs 123 index_t node0;
159 jgs 82 Finley_Mesh* out;
160     char name[50];
161     double time0=Finley_timer();
162     NE0=MAX(1,numElements[0]);
163     NE1=MAX(1,numElements[1]);
164     NE2=MAX(1,numElements[2]);
165     N0=NE0+1;
166     N1=NE1+1;
167     N2=NE2+1;
168 bcumming 751 #ifdef PASO_MPI
169     Paso_MPIInfo *mpi_info = NULL;
170 jgs 82
171 bcumming 751 /* get MPI information */
172     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
173     if (! Finley_noError())
174     return NULL;
175     #endif
176    
177 jgs 153 if (N0<=MIN(N1,N2)) {
178     if (N1 <= N2) {
179     M0=1;
180     M1=N0;
181     M2=N0*N1;
182     } else {
183     M0=1;
184     M2=N0;
185     M1=N0*N2;
186     }
187     } else if (N1<=MIN(N2,N0)) {
188     if (N2 <= N0) {
189     M1=1;
190     M2=N1;
191     M0=N2*N1;
192     } else {
193     M1=1;
194     M0=N1;
195     M2=N1*N0;
196     }
197     } else {
198     if (N0 <= N1) {
199     M2=1;
200     M0=N2;
201     M1=N2*N0;
202     } else {
203     M2=1;
204     M1=N2;
205     M0=N1*N2;
206     }
207     }
208    
209    
210 jgs 82 NFaceElements=0;
211     if (!periodic[0]) {
212     NDOF0=N0;
213     NFaceElements+=2*NE1*NE2;
214     } else {
215     NDOF0=N0-1;
216     }
217     if (!periodic[1]) {
218     NDOF1=N1;
219     NFaceElements+=2*NE0*NE2;
220     } else {
221     NDOF1=N1-1;
222     }
223     if (!periodic[2]) {
224     NDOF2=N2;
225     NFaceElements+=2*NE0*NE1;
226     } else {
227     NDOF2=N2-1;
228     }
229    
230     /* allocate mesh: */
231    
232     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
233 bcumming 751
234 bcumming 730 #ifndef PASO_MPI
235 jgs 82 out=Finley_Mesh_alloc(name,3,order);
236 bcumming 730 #else
237 bcumming 751 out=Finley_Mesh_alloc(name,3,order,mpi_info);
238 bcumming 730 #endif
239 jgs 150 if (! Finley_noError()) return NULL;
240 jgs 82
241 bcumming 751 #ifdef PASO_MPI
242     out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
243     if (useElementsOnFace) {
244     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
245     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
246     } else {
247     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
248     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
249     }
250     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
251     #else
252 jgs 82 out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
253     if (useElementsOnFace) {
254     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
255     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);
256     } else {
257     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);
258     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
259     }
260     out->Points=Finley_ElementFile_alloc(Point1,out->order);
261 bcumming 730 #endif
262 jgs 150 if (! Finley_noError()) {
263 jgs 82 Finley_Mesh_dealloc(out);
264     return NULL;
265     }
266    
267    
268 bcumming 751 /* allocate tables: */
269 jgs 82 Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
270 bcumming 751 #ifdef PASO_MPI
271     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
272 bcumming 730 #endif
273 jgs 82 Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
274     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
275 jgs 150 if (! Finley_noError()) {
276 jgs 82 Finley_Mesh_dealloc(out);
277     return NULL;
278     }
279    
280     /* set nodes: */
281 jgs 153
282     #pragma omp parallel for private(i0,i1,i2,k)
283 jgs 82 for (i2=0;i2<N2;i2++) {
284     for (i1=0;i1<N1;i1++) {
285     for (i0=0;i0<N0;i0++) {
286 jgs 153 k=M0*i0+M1*i1+M2*i2;
287 jgs 82 out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
288     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
289     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
290 jgs 153 out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
291 jgs 82 out->Nodes->Tag[k]=0;
292 jgs 153 out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
293 jgs 82 }
294     }
295     }
296     /* tags for the faces: */
297 jgs 153 /* tags for the faces: */
298 jgs 82 if (!periodic[2]) {
299 jgs 153 for (i1=0;i1<N1;i1++) {
300     for (i0=0;i0<N0;i0++) {
301     out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
302     out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
303     }
304     }
305 jgs 82 }
306     if (!periodic[1]) {
307     for (i2=0;i2<N2;i2++) {
308     for (i0=0;i0<N0;i0++) {
309 jgs 153 out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
310     out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
311 jgs 82 }
312     }
313     }
314     if (!periodic[0]) {
315     for (i2=0;i2<N2;i2++) {
316     for (i1=0;i1<N1;i1++) {
317 jgs 153 out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
318     out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
319 jgs 82 }
320     }
321     }
322 jgs 153
323 jgs 82 /* set the elements: */
324    
325     #pragma omp parallel for private(i0,i1,i2,k,node0)
326     for (i2=0;i2<NE2;i2++) {
327     for (i1=0;i1<NE1;i1++) {
328     for (i0=0;i0<NE0;i0++) {
329     k=i0+NE0*i1+NE0*NE1*i2;
330     node0=i0+i1*N0+N0*N1*i2;
331    
332     out->Elements->Id[k]=k;
333     out->Elements->Tag[k]=0;
334     out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
335    
336     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
337     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
338     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;
339     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;
340     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;
341     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;
342     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;
343     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;
344    
345     }
346     }
347     }
348 jgs 123 out->Elements->minColor=0;
349     out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
350 jgs 82
351     /* face elements: */
352    
353     if (useElementsOnFace) {
354     NUMNODES=8;
355     } else {
356     NUMNODES=4;
357     }
358     totalNECount=NE0*NE1*NE2;
359     faceNECount=0;
360    
361     /* these are the quadrilateral elements on boundary 1 (x3=0): */
362     if (!periodic[2]) {
363     /* ** elements on boundary 100 (x3=0): */
364    
365     #pragma omp parallel for private(i0,i1,k,node0)
366     for (i1=0;i1<NE1;i1++) {
367     for (i0=0;i0<NE0;i0++) {
368     k=i0+NE0*i1+faceNECount;
369     node0=i0+i1*N0;
370    
371     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
372     out->FaceElements->Tag[k]=100;
373     out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
374    
375     if (useElementsOnFace) {
376     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
377     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
378     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
379     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
380     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;
381     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;
382     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
383     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;
384     } else {
385     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
386     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;
387     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;
388     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
389     }
390     }
391     }
392     totalNECount+=NE1*NE0;
393     faceNECount+=NE1*NE0;
394    
395     /* ** elements on boundary 200 (x3=1) */
396    
397     #pragma omp parallel for private(i0,i1,k,node0)
398     for (i1=0;i1<NE1;i1++) {
399     for (i0=0;i0<NE0;i0++) {
400     k=i0+NE0*i1+faceNECount;
401     node0=i0+i1*N0+N0*N1*(NE2-1);
402    
403     out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
404     out->FaceElements->Tag[k]=200;
405     out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
406    
407     if (useElementsOnFace) {
408     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
409     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
410     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
411     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
412     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
413     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
414     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;
415     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;
416     } else {
417     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
418     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;
419     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;
420     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;
421     }
422    
423    
424     }
425     }
426     totalNECount+=NE1*NE0;
427     faceNECount+=NE1*NE0;
428     }
429     if (!periodic[0]) {
430     /* ** elements on boundary 001 (x1=0): */
431    
432     #pragma omp parallel for private(i1,i2,k,node0)
433     for (i2=0;i2<NE2;i2++) {
434     for (i1=0;i1<NE1;i1++) {
435     k=i1+NE1*i2+faceNECount;
436     node0=i1*N0+N0*N1*i2;
437    
438     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
439     out->FaceElements->Tag[k]=1;
440     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
441    
442     if (useElementsOnFace) {
443     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
444     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
445     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
446     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
447     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
448     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;
449     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;
450     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;
451     } else {
452     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
453     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;
454     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;
455     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;
456     }
457     }
458     }
459     totalNECount+=NE1*NE2;
460     faceNECount+=NE1*NE2;
461    
462     /* ** elements on boundary 002 (x1=1): */
463    
464     #pragma omp parallel for private(i1,i2,k,node0)
465     for (i2=0;i2<NE2;i2++) {
466     for (i1=0;i1<NE1;i1++) {
467     k=i1+NE1*i2+faceNECount;
468     node0=(NE0-1)+i1*N0+N0*N1*i2 ;
469    
470     out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
471     out->FaceElements->Tag[k]=2;
472     out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
473    
474     if (useElementsOnFace) {
475     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
476     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
477     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
478     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
479     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
480     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;
481     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;
482     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;
483     } else {
484     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
485     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;
486     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
487     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;
488     }
489     }
490     }
491     totalNECount+=NE1*NE2;
492     faceNECount+=NE1*NE2;
493     }
494     if (!periodic[1]) {
495     /* ** elements on boundary 010 (x2=0): */
496    
497     #pragma omp parallel for private(i0,i2,k,node0)
498     for (i2=0;i2<NE2;i2++) {
499     for (i0=0;i0<NE0;i0++) {
500     k=i0+NE0*i2+faceNECount;
501     node0=i0+N0*N1*i2;
502    
503     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
504     out->FaceElements->Tag[k]=10;
505     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
506    
507     if (useElementsOnFace) {
508     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
509     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
510     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
511     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
512     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;
513     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;
514     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;
515     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;
516     } else {
517     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
518     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
519     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;
520     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;
521     }
522     }
523     }
524     totalNECount+=NE0*NE2;
525     faceNECount+=NE0*NE2;
526    
527     /* ** elements on boundary 020 (x2=1): */
528    
529     #pragma omp parallel for private(i0,i2,k,node0)
530     for (i2=0;i2<NE2;i2++) {
531     for (i0=0;i0<NE0;i0++) {
532     k=i0+NE0*i2+faceNECount;
533     node0=i0+(NE1-1)*N0+N0*N1*i2;
534    
535     out->FaceElements->Tag[k]=20;
536     out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
537     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
538    
539     if (useElementsOnFace) {
540     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
541     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
542     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
543     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
544     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
545     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;
546     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;
547     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
548     } else {
549     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
550     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;
551     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;
552     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;
553     }
554    
555     }
556     }
557     totalNECount+=NE0*NE2;
558     faceNECount+=NE0*NE2;
559     }
560 jgs 123 out->FaceElements->minColor=0;
561     out->FaceElements->maxColor=23;
562 jgs 82
563     /* face elements done: */
564    
565 bcumming 751 #ifdef PASO_MPI
566     /* make sure that the trivial distribution data is correct */
567     out->FaceElements->elementDistribution->numBoundary = 0;
568     out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = faceNECount;
569     out->Elements->elementDistribution->numBoundary = 0;
570     out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal = out->Elements->numElements;
571     out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
572     out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
573    
574     out->Nodes->degreeOfFreedomDistribution->numInternal = out->Nodes->degreeOfFreedomDistribution->numLocal;
575     out->Nodes->degreeOfFreedomDistribution->numBoundary = 0;
576     #endif
577 jgs 82 /* condense the nodes: */
578     Finley_Mesh_resolveNodeIds(out);
579    
580 bcumming 751 #ifdef PASO_MPI
581     /* setup the CommBuffer */
582     Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
583     if ( !Finley_MPI_noError( mpi_info )) {
584     if( Finley_noError() )
585     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
586     Paso_MPIInfo_dealloc( mpi_info );
587     Finley_Mesh_dealloc(out);
588     return NULL;
589     }
590     #endif
591    
592 jgs 82 /* prepare mesh for further calculatuions:*/
593     Finley_Mesh_prepare(out) ;
594    
595 jgs 149 #ifdef Finley_TRACE
596 jgs 82 printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
597 jgs 149 #endif
598 jgs 82
599 jgs 150 if (! Finley_noError()) {
600 jgs 82 Finley_Mesh_dealloc(out);
601     return NULL;
602     }
603     return out;
604     }
605 bcumming 751
606     #ifdef PASO_MPI
607     Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
608     {
609     dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
610     dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
611     bool_t dom_left, dom_right, dom_internal;
612    
613     index_t N0dom;
614     index_t firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1, node2;
615     index_t targetDomain=-1;
616     bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
617     index_t *indexForward=NULL;
618     Finley_Mesh* out;
619    
620     char name[50];
621     Paso_MPIInfo *mpi_info = NULL;
622     double time0=Finley_timer();
623    
624     NE0=MAX(1,numElements[0]);
625     NE1=MAX(1,numElements[1]);
626     NE2=MAX(1,numElements[2]);
627     N0=NE0+1;
628     N1=NE1+1;
629     N2=NE2+1;
630    
631    
632     /* get MPI information */
633     mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
634     if (! Finley_noError())
635     return NULL;
636    
637     /* use the serial version to generate the mesh for the 1-CPU case */
638     if( mpi_info->size==1 )
639     {
640     Paso_MPIInfo_dealloc( mpi_info );
641     out = Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace );
642     //print_mesh_statistics( out );
643     return out;
644     }
645    
646     if( mpi_info->rank==0 )
647     domLeft = TRUE;
648     if( mpi_info->rank==mpi_info->size-1 )
649     domRight = TRUE;
650     if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
651     domInternal = TRUE;
652    
653     /* dimensions of the local subdomain */
654     domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );
655    
656     /* count Degrees of Freedom along each axis */
657     NDOF0 = N0 - periodic[0];
658     NDOF1 = N1 - periodic[1];
659     NDOF2 = N2 - periodic[2];
660    
661     /* count face elements */
662     /* internal face elements */
663     NFaceElements = 0;
664     if( !periodic[0] )
665     NFaceElements += (domLeft+domRight)*NE1*NE2;
666     if( !periodic[1] )
667     NFaceElements += 2*numElementsInternal*NE2;
668     if( !periodic[2] )
669     NFaceElements += 2*numElementsInternal*NE1;
670     /* boundary face elements */
671     /* this is looks nasty, but it beats a bunch of nested if/then/else carry-on */
672     NFaceElements += 2*( 2 - (domLeft + domRight)*(!periodic[0]) )*( (!periodic[1])*NE2 + (!periodic[2])*NE1 );
673    
674    
675     /* allocate mesh: */
676     sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
677    
678     out=Finley_Mesh_alloc(name,3,order,mpi_info);
679     if (! Finley_noError()) return NULL;
680    
681     out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
682     if (useElementsOnFace) {
683     out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
684     out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
685     } else {
686     out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
687     out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
688     }
689     out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
690    
691     if (! Finley_noError()) {
692     Finley_Mesh_dealloc(out);
693     return NULL;
694     }
695    
696    
697     /* allocate tables: */
698     Finley_NodeFile_allocTable(out->Nodes,(numNodesLocal+2-!periodic[0]*(domLeft+domRight))*N1*N2);
699     Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, (DOFExternal[0]+DOFExternal[1])*NDOF1*NDOF2, 0 );
700     Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
701     Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
702     if (! Finley_noError()) {
703     Finley_Mesh_dealloc(out);
704     return NULL;
705     }
706    
707     /* set nodes: */
708     /* INTERNAL/BOUNDARY NODES */
709 bcumming 759 k=0;
710 bcumming 751 #pragma omp parallel for private(i0,i1,i2,k)
711 bcumming 759 for (i2=0;i2<N2;i2++) {
712 bcumming 751 for (i1=0;i1<N1;i1++) {
713     for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {
714     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
715     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
716     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
717     out->Nodes->Id[k]=k;
718     out->Nodes->Tag[k]=0;
719     out->Nodes->degreeOfFreedom[k]=k;
720     }
721     }
722     }
723     if( domLeft && periodic[0] ) {
724     for (i2=0;i2<N2;i2++) {
725     for (i1=0;i1<N1;i1++, k++) {
726     out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
727     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
728     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
729     out->Nodes->Id[k]=k;
730     out->Nodes->Tag[k]=0;
731     out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;
732     }
733     }
734     /* tag the faces for this special case */
735     if( !periodic[1] )
736     {
737     for( i2=0; i2<N2; i2++ ){
738     out->Nodes->Tag[k + (i2-N2)*N1 ] += 10;
739     out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;
740     }
741     }
742     if( !periodic[2] )
743     {
744     for( i1=0; i1<N1; i1++ ){
745     out->Nodes->Tag[k -N1*N2 +i1] += 100;
746     out->Nodes->Tag[k -N1 +i1] += 200;
747     }
748     }
749     }
750     /* tags for the faces: */
751     N0dom = (numNodesLocal-periodicLocal[0]);
752     if (!periodic[2]) {
753     for (i1=0;i1<N1;i1++) {
754     for (i0=0;i0<N0dom;i0++) {
755     out->Nodes->Tag[i0 + N0dom*i1]+=100;
756     out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;
757     }
758     }
759     }
760     if (!periodic[1]) {
761     for (i2=0;i2<N2;i2++) {
762     for (i0=0;i0<N0dom;i0++) {
763     out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;
764     out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;
765     }
766     }
767     }
768     if (!periodic[0] && !domInternal ) {
769     for (i2=0;i2<N2;i2++) {
770     for (i1=0;i1<N1;i1++) {
771     if( domLeft )
772     out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;
773     if( domRight )
774     out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;
775     }
776     }
777     }
778     /* setup the forward communication data for the boundary nodes that we have just defined */
779     /* the case where there are 2 subdomains and periodic[0]=true has to be treated
780     as a special case to because the two domains have two interface boundaries to one-another */
781     indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );
782     if( mpi_info->size>2 || !periodic[0] ){
783     if( domInternal || domRight || periodicLocal[0] )
784     {
785     for( int i=0; i<NDOF2; i++ )
786     for( int j=0; j<NDOF1; j++ )
787     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
788     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
789     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
790     }
791     if( domInternal || domLeft || periodicLocal[1] )
792     {
793     for( int i=0; i<NDOF2; i++ )
794     for( int j=0; j<NDOF1; j++ )
795     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
796     targetDomain = (mpi_info->rank+1) % mpi_info->size;
797     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
798     }
799     }
800     else {
801     for( int i=0; i<NDOF2; i++ )
802     for( int j=0; j<NDOF1; j++ )
803     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
804     targetDomain = (mpi_info->rank+1) % mpi_info->size;
805     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
806    
807     for( int i=0; i<NDOF2; i++ )
808     for( int j=0; j<NDOF1; j++ )
809     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
810     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
811     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
812     }
813    
814     /* EXTERNAL NODES */
815     /* left hand boundary */
816     DOFcount = NDOF1*NDOF2*numDOFLocal;
817     if( (domLeft && periodic[0]) || !domLeft ) {
818     if( (domLeft && periodic[0]) )
819     for (i2=0;i2<N2;i2++) {
820     for (i1=0;i1<N1;i1++, k++) {
821     out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];
822     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
823     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
824     out->Nodes->Id[k]=k;
825     out->Nodes->Tag[k]=0;
826     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
827     }
828     }
829     else
830     for (i2=0;i2<N2;i2++) {
831     for (i1=0;i1<N1;i1++, k++) {
832     out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];
833     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
834     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
835     out->Nodes->Id[k]=k;
836     out->Nodes->Tag[k]=0;
837     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
838     }
839     }
840     DOFcount += NDOF1*NDOF2;
841     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
842     if( !periodic[1] ){
843     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
844     }
845     else {
846     for( int i=0; i<NDOF2; i++ )
847     for( int j=0; j<NDOF1; j++ )
848     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
849     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
850     }
851    
852     /* tag the faces for this special case */
853     if( !periodic[1] )
854     {
855     for( i1=0; i1<N1; i1++ ){
856     out->Nodes->Tag[k -N1*N2 +i1] += 10;
857     out->Nodes->Tag[k -N1 +i1] += 20;
858     }
859     }
860     if( periodic[2] )
861     {
862     for( i2=0; i2<N2; i2++ ){
863     out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
864     out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
865     }
866     }
867     }
868     if( (domRight && periodic[0]) || !domRight )
869     {
870     if( domRight && periodic[0] )
871     for (i2=0;i2<N2;i2++) {
872     for (i1=0;i1<N1;i1++, k++) {
873     out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
874     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
875     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
876     out->Nodes->Id[k]=k;
877     out->Nodes->Tag[k]=0;
878     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
879     }
880     }
881     else
882     for (i2=0;i2<N2;i2++) {
883     for (i1=0;i1<N1;i1++, k++) {
884     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
885     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
886     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
887     out->Nodes->Id[k]=k;
888     out->Nodes->Tag[k]=0;
889     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
890     }
891     }
892     DOFcount += NDOF1*NDOF2;
893    
894     targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;
895     if( !periodic[1] ){
896     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
897     }
898     else {
899     for( int i=0; i<NDOF2; i++ )
900     for( int j=0; j<NDOF1; j++ )
901     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
902     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
903     }
904    
905     /* tag the faces for this special case */
906     if( !periodic[1] )
907     {
908     for( i1=0; i1<N1; i1++ ){
909     out->Nodes->Tag[k -N1*N2 +i1] += 10;
910     out->Nodes->Tag[k -N1 +i1] += 20;
911     }
912     }
913     if( !periodic[2] )
914     {
915     for( i2=0; i2<N2; i2++ ){
916     out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
917     out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
918     }
919     }
920     }
921     out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));
922     out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;
923    
924     TMPMEMFREE( indexForward );
925     /* set the elements: */
926    
927     /* INTERNAL elements */
928     N0dom = (numNodesLocal-periodicLocal[0]);
929     k = 0;
930     #pragma omp parallel for private(i0,i1,i2,k,node0)
931     for (i2=0;i2<NE2;i2++) {
932     for (i1=0;i1<NE1;i1++) {
933     for (i0=0;i0<numElementsInternal;i0++,k++) {
934     node0=i0+i1*N0dom+N0dom*N1*i2;
935    
936     out->Elements->Id[k]=k;
937    
938     out->Elements->Tag[k]=0;
939     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
940    
941     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
942     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
943     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;
944     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;
945     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;
946     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;
947     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;
948     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;
949    
950     }
951     }
952     }
953     out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;
954     out->Elements->elementDistribution->numBoundary = 0;
955    
956     /* BOUNDARY Elements */
957     /* left boundary */
958     if( !domLeft )
959     {
960     for (i2=0;i2<NE2;i2++) {
961     node0 = numNodesLocal*N1*N2 + i2*N1;
962     for (i1=0;i1<NE1;i1++,node0++,k++) {
963     out->Elements->Id[k]=k;
964     out->Elements->Tag[k]=0;
965     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
966    
967     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
968     out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;
969     out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;
970     out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;
971     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;
972     out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;
973     out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;
974     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;
975     }
976     }
977     out->Elements->elementDistribution->numBoundary += NE1*NE2;
978     }
979     /* the left periodic boundary is done a little differently to a left internal boundary */
980     else if( (domLeft && periodic[0]) )
981     {
982     for (i2=0;i2<NE2;i2++) {
983     node0 = numDOFLocal*N1*N2 + i2*N1;
984     node1 = node0 + N1*N2;
985     for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {
986    
987     out->Elements->Id[k]=k;
988     out->Elements->Tag[k]=0;
989     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
990    
991     out->Elements->Nodes[INDEX2(0,k,8)]=node1;
992     out->Elements->Nodes[INDEX2(1,k,8)]=node0;
993     out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
994     out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;
995     out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;
996     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
997     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
998     out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;
999     }
1000     }
1001     out->Elements->elementDistribution->numBoundary += NE1*NE2;
1002     }
1003     /* right boundary */
1004     if( !domRight || (domRight && periodic[0]) ){
1005     for (i2=0;i2<NE2;i2++) {
1006     for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {
1007     node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;
1008     node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;
1009    
1010     out->Elements->Id[k]=k;
1011     out->Elements->Tag[k]=0;
1012     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
1013    
1014     out->Elements->Nodes[INDEX2(0,k,8)]=node1;
1015     out->Elements->Nodes[INDEX2(1,k,8)]=node0;
1016     out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
1017     out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;
1018     out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;
1019     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
1020     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
1021     out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;
1022     }
1023     }
1024     out->Elements->elementDistribution->numBoundary += NE1*NE2;
1025     }
1026    
1027     out->Elements->minColor=0;
1028     out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1029     out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;
1030    
1031     /* face elements: */
1032    
1033     if (useElementsOnFace) {
1034     NUMNODES=8;
1035     } else {
1036     NUMNODES=4;
1037     }
1038     totalNECount=k;
1039     faceNECount=0;
1040     idCount = totalNECount;
1041    
1042     /* Do internal face elements for each boundary face first */
1043     /* these are the quadrilateral elements on boundary 1 (x3=0): */
1044     numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1045     if (!periodic[2]) {
1046     /* elements on boundary 100 (x3=0): */
1047    
1048     #pragma omp parallel for private(i0,i1,k,node0)
1049     for (i1=0;i1<NE1;i1++) {
1050     for (i0=0; i0<numElementsInternal; i0++) {
1051     k=i0+numElementsInternal*i1+faceNECount;
1052     node0=i0+i1*numDOFLocal;
1053    
1054     out->FaceElements->Id[k]=idCount++;
1055     out->FaceElements->Tag[k]=100;
1056     out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);
1057    
1058     if (useElementsOnFace) {
1059     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1060     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1061     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1062     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1063     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1064     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1065     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1066     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1067     } else {
1068     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1069     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1070     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1071     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1072     }
1073     }
1074     }
1075     totalNECount+=NE1*numElementsInternal;
1076     faceNECount+=NE1*numElementsInternal;
1077    
1078     /* ** elements on boundary 200 (x3=1) */
1079    
1080     #pragma omp parallel for private(i0,i1,k,node0)
1081     for (i1=0;i1<NE1;i1++) {
1082     for (i0=0;i0<numElementsInternal;i0++) {
1083     k=i0+numElementsInternal*i1+faceNECount;
1084     node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);
1085    
1086     out->FaceElements->Id[k]=idCount++;
1087     out->FaceElements->Tag[k]=200;
1088     out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;
1089    
1090     if (useElementsOnFace) {
1091     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1092     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1093     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1094     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1095     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1096     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1097     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;
1098     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1099     } else {
1100     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1101     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1102     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1103     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1104     }
1105     }
1106     }
1107     totalNECount+=NE1*numElementsInternal;
1108     faceNECount+=NE1*numElementsInternal;
1109     }
1110     if (!periodic[0] && !domInternal) {
1111     /* ** elements on boundary 001 (x1=0): */
1112     if( domLeft ){
1113     #pragma omp parallel for private(i1,i2,k,node0)
1114     for (i2=0;i2<NE2;i2++) {
1115     for (i1=0;i1<NE1;i1++) {
1116     k=i1+NE1*i2+faceNECount;
1117     node0=i1*numDOFLocal+numDOFLocal*N1*i2;
1118    
1119     out->FaceElements->Id[k]=idCount++;
1120     out->FaceElements->Tag[k]=1;
1121     out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;
1122    
1123     if (useElementsOnFace) {
1124     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1125     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1126     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1127     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1128     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
1129     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1130     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1131     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;
1132     } else {
1133     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1134     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1135     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1136     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1137     }
1138     }
1139     }
1140     totalNECount+=NE1*NE2;
1141     faceNECount+=NE1*NE2;
1142     }
1143     /* ** elements on boundary 002 (x1=1): */
1144     if( domRight ) {
1145     #pragma omp parallel for private(i1,i2,k,node0)
1146     for (i2=0;i2<NE2;i2++) {
1147     for (i1=0;i1<NE1;i1++) {
1148     k=i1+NE1*i2+faceNECount;
1149     node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;
1150    
1151     out->FaceElements->Id[k]=idCount++;
1152     out->FaceElements->Tag[k]=2;
1153     out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;
1154    
1155     if (useElementsOnFace) {
1156     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1157     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1158     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1159     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1160     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1161     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;
1162     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1163     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;
1164     } else {
1165     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1166     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1167     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1168     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1169     }
1170     }
1171     }
1172     totalNECount+=NE1*NE2;
1173     faceNECount+=NE1*NE2;
1174     }
1175     }
1176     if (!periodic[1]) {
1177     /* ** elements on boundary 010 (x2=0): */
1178    
1179     #pragma omp parallel for private(i0,i2,k,node0)
1180     for (i2=0;i2<NE2;i2++) {
1181     for (i0=0;i0<numElementsInternal;i0++) {
1182     k=i0+numElementsInternal*i2+faceNECount;
1183     node0=i0+numDOFLocal*N1*i2;
1184    
1185     out->FaceElements->Id[k]=idCount++;
1186     out->FaceElements->Tag[k]=10;
1187     out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;
1188    
1189     if (useElementsOnFace) {
1190     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1191     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1192     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1193     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1194     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1195     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;
1196     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;
1197     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1198     } else {
1199     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1200     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1201     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1202     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1203     }
1204     }
1205     }
1206     totalNECount+=numElementsInternal*NE2;
1207     faceNECount+=numElementsInternal*NE2;
1208    
1209     /* ** elements on boundary 020 (x2=1): */
1210    
1211     #pragma omp parallel for private(i0,i2,k,node0)
1212     for (i2=0;i2<NE2;i2++) {
1213     for (i0=0;i0<numElementsInternal;i0++) {
1214     k=i0+numElementsInternal*i2+faceNECount;
1215     node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;
1216    
1217     out->FaceElements->Tag[k]=20;
1218     out->FaceElements->Id[k]=idCount++;
1219     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1220    
1221     if (useElementsOnFace) {
1222     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1223     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1224     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1225     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1226     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1227     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;
1228     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1229     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
1230     } else {
1231     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1232     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1233     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1234     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1235     }
1236    
1237     }
1238     }
1239     totalNECount+=numElementsInternal*NE2;
1240     faceNECount+=numElementsInternal*NE2;
1241     }
1242     out->FaceElements->elementDistribution->numInternal = faceNECount;
1243    
1244     /* now do the boundary face elements */
1245     /* LHS */
1246     if( !domLeft )
1247     {
1248     if( !periodic[2] ) {
1249    
1250     /* x3=0 */
1251     for( i1=0; i1<NE1; i1++ )
1252     {
1253     k = i1+faceNECount;
1254     node0 = i1*numNodesLocal;
1255     node1 = numNodesLocal*N1*N2 + i1;
1256    
1257     out->FaceElements->Tag[k]=200;
1258     out->FaceElements->Id[k]=idCount++;
1259     out->FaceElements->Color[k]=0;
1260    
1261     if( useElementsOnFace )
1262     {
1263     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1264     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1265     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1266     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1267     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1268     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1269     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1270     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;
1271     } else {
1272     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1273     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1274     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1275     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1276     }
1277     }
1278     faceNECount += NE1;
1279     totalNECount += NE1;
1280    
1281     /* x3=1 */
1282     for( i1=0; i1<NE1; i1++ )
1283     {
1284     k = i1+faceNECount;
1285     node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;
1286     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1287    
1288     out->FaceElements->Tag[k]=200;
1289     out->FaceElements->Id[k]=idCount++;
1290     out->FaceElements->Color[k]=0;
1291    
1292     if( useElementsOnFace )
1293     {
1294     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1295     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1296     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1297     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1298     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1299     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1300     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;
1301     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1302     } else {
1303     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1304     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1305     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1306     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1307     }
1308     }
1309     faceNECount += NE1;
1310     totalNECount += NE1;
1311     }
1312    
1313     if( !periodic[1] ) {
1314     /* x2=0 */
1315     for( i2=0; i2<NE2; i2++ )
1316     {
1317     k = i2+faceNECount;
1318     node0 = i2*numNodesLocal*N1;
1319     node1 = numNodesLocal*N1*N2 + i2*N1;
1320    
1321     out->FaceElements->Tag[k]=20;
1322     out->FaceElements->Id[k]=idCount++;
1323     out->FaceElements->Color[k]=0;
1324    
1325     if( useElementsOnFace )
1326     {
1327     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1328     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1329     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1330     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1331     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1332     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;
1333     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1334     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1335     } else {
1336     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1337     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1338     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1339     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1340     }
1341     }
1342     faceNECount += NE2;
1343     totalNECount += NE2;
1344    
1345     /* x2=1 */
1346     for( i2=0; i2<NE2; i2++ )
1347     {
1348     k = i2+faceNECount;
1349     node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);
1350     node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);
1351    
1352     out->FaceElements->Tag[k]=20;
1353     out->FaceElements->Id[k]=idCount++;
1354     out->FaceElements->Color[k]=0;
1355    
1356     if( useElementsOnFace )
1357     {
1358     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1359     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1360     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;
1361     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1362     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1363     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1364     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1365     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1366     } else {
1367     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1368     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1369     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1370     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1371     }
1372     }
1373     faceNECount += NE2;
1374     totalNECount += NE2;
1375     }
1376     }
1377    
1378     /* RHS */
1379     if( !domRight || periodicLocal[1] )
1380     {
1381     /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */
1382     if( domLeft && periodic[0] ){
1383     if( !periodic[2] ) {
1384    
1385     /* x3=0 */
1386     for( i1=0; i1<NE1; i1++ )
1387     {
1388     k = i1+faceNECount;
1389     node0 = numDOFLocal*N1*N2 + i1;
1390     node1 = numNodesLocal*N1*N2 + i1;
1391    
1392     out->FaceElements->Tag[k]=200;
1393     out->FaceElements->Id[k]=idCount++;
1394     out->FaceElements->Color[k]=0;
1395    
1396     if( useElementsOnFace )
1397     {
1398     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1399     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1400     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1401     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1402     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1403     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1404     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1405     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;
1406     } else {
1407     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1408     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1409     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1410     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1411     }
1412     }
1413     faceNECount += NE1;
1414     totalNECount += NE1;
1415    
1416     /* x3=1 */
1417     for( i1=0; i1<NE1; i1++ )
1418     {
1419     k = i1+faceNECount;
1420     node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;
1421     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1422    
1423     out->FaceElements->Tag[k]=200;
1424     out->FaceElements->Id[k]=idCount++;
1425     out->FaceElements->Color[k]=0;
1426    
1427     if( useElementsOnFace )
1428     {
1429     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1430     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1431     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
1432     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1433     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1434     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1435     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1436     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1437     } else {
1438     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1439     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1440     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1441     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1442     }
1443     }
1444     faceNECount += NE1;
1445     totalNECount += NE1;
1446     }
1447    
1448     if( !periodic[1] ) {
1449     /* x2=0 */
1450     for( i2=0; i2<NE2; i2++ )
1451     {
1452     k = i2+faceNECount;
1453     node0 = numDOFLocal*N1*N2 + i2*N1;
1454     node1 = numNodesLocal*N1*N2 + i2*N1;
1455    
1456     out->FaceElements->Tag[k]=20;
1457     out->FaceElements->Id[k]=idCount++;
1458     out->FaceElements->Color[k]=0;
1459    
1460     if( useElementsOnFace )
1461     {
1462     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1463     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1464     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1465     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1466     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1467     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1468     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1469     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1470     } else {
1471     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1472     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1473     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1474     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1475     }
1476     }
1477     faceNECount += NE2;
1478     totalNECount += NE2;
1479    
1480     /* x2=1 */
1481     for( i2=0; i2<NE2; i2++ )
1482     {
1483     k = i2+faceNECount;
1484     node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;
1485     node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;
1486    
1487     out->FaceElements->Tag[k]=20;
1488     out->FaceElements->Id[k]=idCount++;
1489     out->FaceElements->Color[k]=0;
1490    
1491     if( useElementsOnFace )
1492     {
1493     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1494     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1495     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;
1496     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1497     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1498     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1499     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1500     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1501     } else {
1502     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1503     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1504     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1505     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1506     }
1507     }
1508     faceNECount += NE2;
1509     totalNECount += NE2;
1510     }
1511    
1512     }
1513     if( !periodic[2] ) {
1514     /* x3=0 */
1515     for( i1=0; i1<NE1; i1++ )
1516     {
1517     k = i1+faceNECount;
1518     node0 = numDOFLocal*(i1+1) - 1;
1519     node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;
1520    
1521     out->FaceElements->Tag[k]=200;
1522     out->FaceElements->Id[k]=idCount++;
1523     out->FaceElements->Color[k]=0;
1524    
1525     if( useElementsOnFace )
1526     {
1527     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1528     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1529     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1530     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1531     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1532     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1533     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1534     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;
1535     } else {
1536     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1537     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1538     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1539     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1540     }
1541     }
1542     faceNECount += NE1;
1543     totalNECount += NE1;
1544    
1545     /* x3=1 */
1546     for( i1=0; i1<NE1; i1++ )
1547     {
1548     k = i1+faceNECount;
1549     node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;
1550     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;
1551    
1552     out->FaceElements->Tag[k]=200;
1553     out->FaceElements->Id[k]=idCount++;
1554     out->FaceElements->Color[k]=0;
1555    
1556     if( useElementsOnFace )
1557     {
1558     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1559     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1560     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;
1561     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;
1562     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1563     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1564     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1565     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1566     } else {
1567     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1568     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1569     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1570     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1571     }
1572     }
1573     faceNECount += NE1;
1574     totalNECount += NE1;
1575     }
1576     if( !periodic[1] ) {
1577     /* x2=0 */
1578     for( i2=0; i2<NE2; i2++ )
1579     {
1580     k = i2+faceNECount;
1581     node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;
1582     node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1583    
1584     out->FaceElements->Tag[k]=20;
1585     out->FaceElements->Id[k]=idCount++;
1586     out->FaceElements->Color[k]=0;
1587    
1588     if( useElementsOnFace )
1589     {
1590     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1591     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1592     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1593     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1594     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1595     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;
1596     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1597     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1598     } else {
1599     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1600     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1601     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1602     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1603     }
1604     }
1605     faceNECount += NE2;
1606     totalNECount += NE2;
1607    
1608     /* x2=1 */
1609     for( i2=0; i2<NE2; i2++ )
1610     {
1611     k = i2+faceNECount;
1612     node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;
1613     node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1614    
1615     out->FaceElements->Tag[k]=20;
1616     out->FaceElements->Id[k]=idCount++;
1617     out->FaceElements->Color[k]=0;
1618    
1619     if( useElementsOnFace ){
1620     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1621     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;
1622     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;
1623     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;
1624     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1625     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1626     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1627     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1628     } else {
1629     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1630     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1631     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1632     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1633     }
1634     }
1635     faceNECount += NE2;
1636     totalNECount += NE2;
1637     }
1638     }
1639     out->FaceElements->minColor=0;
1640     out->FaceElements->maxColor=0;//23;
1641    
1642     out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;
1643     out->FaceElements->elementDistribution->numLocal = faceNECount;
1644    
1645    
1646     /* setup distribution info for other elements */
1647     out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
1648     out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
1649    
1650     /* condense the nodes: */
1651     Finley_Mesh_resolveNodeIds( out );
1652    
1653     /* setup the CommBuffer */
1654     Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1655     if ( !Finley_MPI_noError( mpi_info )) {
1656     if( Finley_noError() )
1657     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1658     Paso_MPIInfo_dealloc( mpi_info );
1659     Finley_Mesh_dealloc(out);
1660     return NULL;
1661     }
1662    
1663     Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1664    
1665     /* prepare mesh for further calculatuions:*/
1666     Finley_Mesh_prepare(out) ;
1667    
1668     // print_mesh_statistics( out );
1669    
1670     #ifdef Finley_TRACE
1671     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1672     #endif
1673    
1674     if (! Finley_noError()) {
1675     Finley_Mesh_dealloc(out);
1676     return NULL;
1677     }
1678    
1679     return out;
1680     }
1681     #endif
1682    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26