/[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 751 - (hide annotations)
Mon Jun 26 01:46:34 2006 UTC (13 years ago) by bcumming
File MIME type: text/plain
File size: 64654 byte(s)
Changes relating to the MPI version of escript
The standard OpenMP version of escript is unchanged

- updated data types (Finley_Mesh, Finley_NodeFile, etc) to store meshes
  over multiple MPI processes.
- added CommBuffer code in Paso for communication of Data associated
  with distributed meshes
- updates in Finley and Escript to support distributed data and operations
  on distributed data (such as interpolation).
- construction of RHS in MPI, so that simple explicit schemes (such as
  /docs/examples/wave.py without IO and the Locator) can run in MPI.
- updated mesh generation for first order line, rectangle and brick
  meshes and second order line meshes in MPI.        
- small changes to trunk/SConstruct and trunk/scons/ess_options.py to
  build the MPI version, these changes are turned off by default.

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     #pragma omp parallel for private(i0,i1,i2,k)
710     for (k=0,i2=0;i2<N2;i2++) {
711     for (i1=0;i1<N1;i1++) {
712     for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {
713     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
714     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
715     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
716     out->Nodes->Id[k]=k;
717     out->Nodes->Tag[k]=0;
718     out->Nodes->degreeOfFreedom[k]=k;
719     }
720     }
721     }
722     if( domLeft && periodic[0] ) {
723     for (i2=0;i2<N2;i2++) {
724     for (i1=0;i1<N1;i1++, k++) {
725     out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
726     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
727     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
728     out->Nodes->Id[k]=k;
729     out->Nodes->Tag[k]=0;
730     out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;
731     }
732     }
733     /* tag the faces for this special case */
734     if( !periodic[1] )
735     {
736     for( i2=0; i2<N2; i2++ ){
737     out->Nodes->Tag[k + (i2-N2)*N1 ] += 10;
738     out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;
739     }
740     }
741     if( !periodic[2] )
742     {
743     for( i1=0; i1<N1; i1++ ){
744     out->Nodes->Tag[k -N1*N2 +i1] += 100;
745     out->Nodes->Tag[k -N1 +i1] += 200;
746     }
747     }
748     }
749     /* tags for the faces: */
750     N0dom = (numNodesLocal-periodicLocal[0]);
751     if (!periodic[2]) {
752     for (i1=0;i1<N1;i1++) {
753     for (i0=0;i0<N0dom;i0++) {
754     out->Nodes->Tag[i0 + N0dom*i1]+=100;
755     out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;
756     }
757     }
758     }
759     if (!periodic[1]) {
760     for (i2=0;i2<N2;i2++) {
761     for (i0=0;i0<N0dom;i0++) {
762     out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;
763     out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;
764     }
765     }
766     }
767     if (!periodic[0] && !domInternal ) {
768     for (i2=0;i2<N2;i2++) {
769     for (i1=0;i1<N1;i1++) {
770     if( domLeft )
771     out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;
772     if( domRight )
773     out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;
774     }
775     }
776     }
777     /* setup the forward communication data for the boundary nodes that we have just defined */
778     /* the case where there are 2 subdomains and periodic[0]=true has to be treated
779     as a special case to because the two domains have two interface boundaries to one-another */
780     indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );
781     if( mpi_info->size>2 || !periodic[0] ){
782     if( domInternal || domRight || periodicLocal[0] )
783     {
784     for( int i=0; i<NDOF2; i++ )
785     for( int j=0; j<NDOF1; j++ )
786     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
787     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
788     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
789     }
790     if( domInternal || domLeft || periodicLocal[1] )
791     {
792     for( int i=0; i<NDOF2; i++ )
793     for( int j=0; j<NDOF1; j++ )
794     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
795     targetDomain = (mpi_info->rank+1) % mpi_info->size;
796     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
797     }
798     }
799     else {
800     for( int i=0; i<NDOF2; i++ )
801     for( int j=0; j<NDOF1; j++ )
802     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
803     targetDomain = (mpi_info->rank+1) % mpi_info->size;
804     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
805    
806     for( int i=0; i<NDOF2; i++ )
807     for( int j=0; j<NDOF1; j++ )
808     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
809     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
810     Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
811     }
812    
813     /* EXTERNAL NODES */
814     /* left hand boundary */
815     DOFcount = NDOF1*NDOF2*numDOFLocal;
816     if( (domLeft && periodic[0]) || !domLeft ) {
817     if( (domLeft && periodic[0]) )
818     for (i2=0;i2<N2;i2++) {
819     for (i1=0;i1<N1;i1++, k++) {
820     out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];
821     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
822     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
823     out->Nodes->Id[k]=k;
824     out->Nodes->Tag[k]=0;
825     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
826     }
827     }
828     else
829     for (i2=0;i2<N2;i2++) {
830     for (i1=0;i1<N1;i1++, k++) {
831     out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];
832     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
833     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
834     out->Nodes->Id[k]=k;
835     out->Nodes->Tag[k]=0;
836     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
837     }
838     }
839     DOFcount += NDOF1*NDOF2;
840     targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
841     if( !periodic[1] ){
842     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
843     }
844     else {
845     for( int i=0; i<NDOF2; i++ )
846     for( int j=0; j<NDOF1; j++ )
847     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
848     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
849     }
850    
851     /* tag the faces for this special case */
852     if( !periodic[1] )
853     {
854     for( i1=0; i1<N1; i1++ ){
855     out->Nodes->Tag[k -N1*N2 +i1] += 10;
856     out->Nodes->Tag[k -N1 +i1] += 20;
857     }
858     }
859     if( periodic[2] )
860     {
861     for( i2=0; i2<N2; i2++ ){
862     out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
863     out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
864     }
865     }
866     }
867     if( (domRight && periodic[0]) || !domRight )
868     {
869     if( domRight && periodic[0] )
870     for (i2=0;i2<N2;i2++) {
871     for (i1=0;i1<N1;i1++, k++) {
872     out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
873     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
874     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
875     out->Nodes->Id[k]=k;
876     out->Nodes->Tag[k]=0;
877     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
878     }
879     }
880     else
881     for (i2=0;i2<N2;i2++) {
882     for (i1=0;i1<N1;i1++, k++) {
883     out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
884     out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
885     out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
886     out->Nodes->Id[k]=k;
887     out->Nodes->Tag[k]=0;
888     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
889     }
890     }
891     DOFcount += NDOF1*NDOF2;
892    
893     targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;
894     if( !periodic[1] ){
895     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
896     }
897     else {
898     for( int i=0; i<NDOF2; i++ )
899     for( int j=0; j<NDOF1; j++ )
900     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
901     Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
902     }
903    
904     /* tag the faces for this special case */
905     if( !periodic[1] )
906     {
907     for( i1=0; i1<N1; i1++ ){
908     out->Nodes->Tag[k -N1*N2 +i1] += 10;
909     out->Nodes->Tag[k -N1 +i1] += 20;
910     }
911     }
912     if( !periodic[2] )
913     {
914     for( i2=0; i2<N2; i2++ ){
915     out->Nodes->Tag[k +(i2-N2)*N1 ] += 100;
916     out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;
917     }
918     }
919     }
920     out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));
921     out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;
922    
923     TMPMEMFREE( indexForward );
924     /* set the elements: */
925    
926     /* INTERNAL elements */
927     N0dom = (numNodesLocal-periodicLocal[0]);
928     k = 0;
929     #pragma omp parallel for private(i0,i1,i2,k,node0)
930     for (i2=0;i2<NE2;i2++) {
931     for (i1=0;i1<NE1;i1++) {
932     for (i0=0;i0<numElementsInternal;i0++,k++) {
933     node0=i0+i1*N0dom+N0dom*N1*i2;
934    
935     out->Elements->Id[k]=k;
936    
937     out->Elements->Tag[k]=0;
938     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
939    
940     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
941     out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
942     out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;
943     out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;
944     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;
945     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;
946     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;
947     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;
948    
949     }
950     }
951     }
952     out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;
953     out->Elements->elementDistribution->numBoundary = 0;
954    
955     /* BOUNDARY Elements */
956     /* left boundary */
957     if( !domLeft )
958     {
959     for (i2=0;i2<NE2;i2++) {
960     node0 = numNodesLocal*N1*N2 + i2*N1;
961     for (i1=0;i1<NE1;i1++,node0++,k++) {
962     out->Elements->Id[k]=k;
963     out->Elements->Tag[k]=0;
964     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
965    
966     out->Elements->Nodes[INDEX2(0,k,8)]=node0;
967     out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;
968     out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;
969     out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;
970     out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;
971     out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;
972     out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;
973     out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;
974     }
975     }
976     out->Elements->elementDistribution->numBoundary += NE1*NE2;
977     }
978     /* the left periodic boundary is done a little differently to a left internal boundary */
979     else if( (domLeft && periodic[0]) )
980     {
981     for (i2=0;i2<NE2;i2++) {
982     node0 = numDOFLocal*N1*N2 + i2*N1;
983     node1 = node0 + N1*N2;
984     for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {
985    
986     out->Elements->Id[k]=k;
987     out->Elements->Tag[k]=0;
988     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
989    
990     out->Elements->Nodes[INDEX2(0,k,8)]=node1;
991     out->Elements->Nodes[INDEX2(1,k,8)]=node0;
992     out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
993     out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;
994     out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;
995     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
996     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
997     out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;
998     }
999     }
1000     out->Elements->elementDistribution->numBoundary += NE1*NE2;
1001     }
1002     /* right boundary */
1003     if( !domRight || (domRight && periodic[0]) ){
1004     for (i2=0;i2<NE2;i2++) {
1005     for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {
1006     node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;
1007     node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;
1008    
1009     out->Elements->Id[k]=k;
1010     out->Elements->Tag[k]=0;
1011     out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
1012    
1013     out->Elements->Nodes[INDEX2(0,k,8)]=node1;
1014     out->Elements->Nodes[INDEX2(1,k,8)]=node0;
1015     out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
1016     out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;
1017     out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;
1018     out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
1019     out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
1020     out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;
1021     }
1022     }
1023     out->Elements->elementDistribution->numBoundary += NE1*NE2;
1024     }
1025    
1026     out->Elements->minColor=0;
1027     out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1028     out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;
1029    
1030     /* face elements: */
1031    
1032     if (useElementsOnFace) {
1033     NUMNODES=8;
1034     } else {
1035     NUMNODES=4;
1036     }
1037     totalNECount=k;
1038     faceNECount=0;
1039     idCount = totalNECount;
1040    
1041     /* Do internal face elements for each boundary face first */
1042     /* these are the quadrilateral elements on boundary 1 (x3=0): */
1043     numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1044     if (!periodic[2]) {
1045     /* elements on boundary 100 (x3=0): */
1046    
1047     #pragma omp parallel for private(i0,i1,k,node0)
1048     for (i1=0;i1<NE1;i1++) {
1049     for (i0=0; i0<numElementsInternal; i0++) {
1050     k=i0+numElementsInternal*i1+faceNECount;
1051     node0=i0+i1*numDOFLocal;
1052    
1053     out->FaceElements->Id[k]=idCount++;
1054     out->FaceElements->Tag[k]=100;
1055     out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);
1056    
1057     if (useElementsOnFace) {
1058     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1059     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1060     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1061     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1062     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1063     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1064     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1065     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1066     } else {
1067     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1068     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1069     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1070     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1071     }
1072     }
1073     }
1074     totalNECount+=NE1*numElementsInternal;
1075     faceNECount+=NE1*numElementsInternal;
1076    
1077     /* ** elements on boundary 200 (x3=1) */
1078    
1079     #pragma omp parallel for private(i0,i1,k,node0)
1080     for (i1=0;i1<NE1;i1++) {
1081     for (i0=0;i0<numElementsInternal;i0++) {
1082     k=i0+numElementsInternal*i1+faceNECount;
1083     node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);
1084    
1085     out->FaceElements->Id[k]=idCount++;
1086     out->FaceElements->Tag[k]=200;
1087     out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;
1088    
1089     if (useElementsOnFace) {
1090     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1091     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1092     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1093     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1094     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1095     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1096     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;
1097     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1098     } else {
1099     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1100     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1101     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1102     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1103     }
1104     }
1105     }
1106     totalNECount+=NE1*numElementsInternal;
1107     faceNECount+=NE1*numElementsInternal;
1108     }
1109     if (!periodic[0] && !domInternal) {
1110     /* ** elements on boundary 001 (x1=0): */
1111     if( domLeft ){
1112     #pragma omp parallel for private(i1,i2,k,node0)
1113     for (i2=0;i2<NE2;i2++) {
1114     for (i1=0;i1<NE1;i1++) {
1115     k=i1+NE1*i2+faceNECount;
1116     node0=i1*numDOFLocal+numDOFLocal*N1*i2;
1117    
1118     out->FaceElements->Id[k]=idCount++;
1119     out->FaceElements->Tag[k]=1;
1120     out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;
1121    
1122     if (useElementsOnFace) {
1123     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1124     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1125     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1126     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1127     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
1128     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1129     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1130     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;
1131     } else {
1132     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1133     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1134     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1135     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1136     }
1137     }
1138     }
1139     totalNECount+=NE1*NE2;
1140     faceNECount+=NE1*NE2;
1141     }
1142     /* ** elements on boundary 002 (x1=1): */
1143     if( domRight ) {
1144     #pragma omp parallel for private(i1,i2,k,node0)
1145     for (i2=0;i2<NE2;i2++) {
1146     for (i1=0;i1<NE1;i1++) {
1147     k=i1+NE1*i2+faceNECount;
1148     node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;
1149    
1150     out->FaceElements->Id[k]=idCount++;
1151     out->FaceElements->Tag[k]=2;
1152     out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;
1153    
1154     if (useElementsOnFace) {
1155     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1156     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1157     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1158     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1159     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1160     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;
1161     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1162     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;
1163     } else {
1164     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1165     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1166     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1167     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1168     }
1169     }
1170     }
1171     totalNECount+=NE1*NE2;
1172     faceNECount+=NE1*NE2;
1173     }
1174     }
1175     if (!periodic[1]) {
1176     /* ** elements on boundary 010 (x2=0): */
1177    
1178     #pragma omp parallel for private(i0,i2,k,node0)
1179     for (i2=0;i2<NE2;i2++) {
1180     for (i0=0;i0<numElementsInternal;i0++) {
1181     k=i0+numElementsInternal*i2+faceNECount;
1182     node0=i0+numDOFLocal*N1*i2;
1183    
1184     out->FaceElements->Id[k]=idCount++;
1185     out->FaceElements->Tag[k]=10;
1186     out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;
1187    
1188     if (useElementsOnFace) {
1189     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1190     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1191     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1192     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1193     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1194     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;
1195     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;
1196     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1197     } else {
1198     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1199     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1200     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1201     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1202     }
1203     }
1204     }
1205     totalNECount+=numElementsInternal*NE2;
1206     faceNECount+=numElementsInternal*NE2;
1207    
1208     /* ** elements on boundary 020 (x2=1): */
1209    
1210     #pragma omp parallel for private(i0,i2,k,node0)
1211     for (i2=0;i2<NE2;i2++) {
1212     for (i0=0;i0<numElementsInternal;i0++) {
1213     k=i0+numElementsInternal*i2+faceNECount;
1214     node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;
1215    
1216     out->FaceElements->Tag[k]=20;
1217     out->FaceElements->Id[k]=idCount++;
1218     out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1219    
1220     if (useElementsOnFace) {
1221     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1222     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1223     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1224     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1225     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1226     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;
1227     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1228     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
1229     } else {
1230     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1231     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1232     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1233     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1234     }
1235    
1236     }
1237     }
1238     totalNECount+=numElementsInternal*NE2;
1239     faceNECount+=numElementsInternal*NE2;
1240     }
1241     out->FaceElements->elementDistribution->numInternal = faceNECount;
1242    
1243     /* now do the boundary face elements */
1244     /* LHS */
1245     if( !domLeft )
1246     {
1247     if( !periodic[2] ) {
1248    
1249     /* x3=0 */
1250     for( i1=0; i1<NE1; i1++ )
1251     {
1252     k = i1+faceNECount;
1253     node0 = i1*numNodesLocal;
1254     node1 = numNodesLocal*N1*N2 + i1;
1255    
1256     out->FaceElements->Tag[k]=200;
1257     out->FaceElements->Id[k]=idCount++;
1258     out->FaceElements->Color[k]=0;
1259    
1260     if( useElementsOnFace )
1261     {
1262     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1263     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1264     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1265     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1266     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1267     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1268     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1269     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;
1270     } else {
1271     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1272     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1273     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1274     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1275     }
1276     }
1277     faceNECount += NE1;
1278     totalNECount += NE1;
1279    
1280     /* x3=1 */
1281     for( i1=0; i1<NE1; i1++ )
1282     {
1283     k = i1+faceNECount;
1284     node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;
1285     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1286    
1287     out->FaceElements->Tag[k]=200;
1288     out->FaceElements->Id[k]=idCount++;
1289     out->FaceElements->Color[k]=0;
1290    
1291     if( useElementsOnFace )
1292     {
1293     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1294     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1295     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1296     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1297     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1298     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1299     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;
1300     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1301     } else {
1302     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1303     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1304     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1305     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1306     }
1307     }
1308     faceNECount += NE1;
1309     totalNECount += NE1;
1310     }
1311    
1312     if( !periodic[1] ) {
1313     /* x2=0 */
1314     for( i2=0; i2<NE2; i2++ )
1315     {
1316     k = i2+faceNECount;
1317     node0 = i2*numNodesLocal*N1;
1318     node1 = numNodesLocal*N1*N2 + i2*N1;
1319    
1320     out->FaceElements->Tag[k]=20;
1321     out->FaceElements->Id[k]=idCount++;
1322     out->FaceElements->Color[k]=0;
1323    
1324     if( useElementsOnFace )
1325     {
1326     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1327     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1328     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1329     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1330     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1331     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;
1332     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1333     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1334     } else {
1335     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1336     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1337     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1338     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1339     }
1340     }
1341     faceNECount += NE2;
1342     totalNECount += NE2;
1343    
1344     /* x2=1 */
1345     for( i2=0; i2<NE2; i2++ )
1346     {
1347     k = i2+faceNECount;
1348     node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);
1349     node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);
1350    
1351     out->FaceElements->Tag[k]=20;
1352     out->FaceElements->Id[k]=idCount++;
1353     out->FaceElements->Color[k]=0;
1354    
1355     if( useElementsOnFace )
1356     {
1357     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1358     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1359     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;
1360     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1361     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1362     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1363     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1364     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1365     } else {
1366     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1367     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1368     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1369     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1370     }
1371     }
1372     faceNECount += NE2;
1373     totalNECount += NE2;
1374     }
1375     }
1376    
1377     /* RHS */
1378     if( !domRight || periodicLocal[1] )
1379     {
1380     /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */
1381     if( domLeft && periodic[0] ){
1382     if( !periodic[2] ) {
1383    
1384     /* x3=0 */
1385     for( i1=0; i1<NE1; i1++ )
1386     {
1387     k = i1+faceNECount;
1388     node0 = numDOFLocal*N1*N2 + i1;
1389     node1 = numNodesLocal*N1*N2 + i1;
1390    
1391     out->FaceElements->Tag[k]=200;
1392     out->FaceElements->Id[k]=idCount++;
1393     out->FaceElements->Color[k]=0;
1394    
1395     if( useElementsOnFace )
1396     {
1397     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1398     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1399     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1400     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1401     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1402     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1403     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1404     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;
1405     } else {
1406     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1407     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1408     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1409     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1410     }
1411     }
1412     faceNECount += NE1;
1413     totalNECount += NE1;
1414    
1415     /* x3=1 */
1416     for( i1=0; i1<NE1; i1++ )
1417     {
1418     k = i1+faceNECount;
1419     node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;
1420     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1421    
1422     out->FaceElements->Tag[k]=200;
1423     out->FaceElements->Id[k]=idCount++;
1424     out->FaceElements->Color[k]=0;
1425    
1426     if( useElementsOnFace )
1427     {
1428     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1429     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1430     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
1431     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1432     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1433     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1434     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1435     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1436     } else {
1437     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1438     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1439     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1440     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1441     }
1442     }
1443     faceNECount += NE1;
1444     totalNECount += NE1;
1445     }
1446    
1447     if( !periodic[1] ) {
1448     /* x2=0 */
1449     for( i2=0; i2<NE2; i2++ )
1450     {
1451     k = i2+faceNECount;
1452     node0 = numDOFLocal*N1*N2 + i2*N1;
1453     node1 = numNodesLocal*N1*N2 + i2*N1;
1454    
1455     out->FaceElements->Tag[k]=20;
1456     out->FaceElements->Id[k]=idCount++;
1457     out->FaceElements->Color[k]=0;
1458    
1459     if( useElementsOnFace )
1460     {
1461     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1462     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1463     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1464     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1465     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1466     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1467     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1468     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1469     } else {
1470     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1471     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1472     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1473     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1474     }
1475     }
1476     faceNECount += NE2;
1477     totalNECount += NE2;
1478    
1479     /* x2=1 */
1480     for( i2=0; i2<NE2; i2++ )
1481     {
1482     k = i2+faceNECount;
1483     node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;
1484     node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;
1485    
1486     out->FaceElements->Tag[k]=20;
1487     out->FaceElements->Id[k]=idCount++;
1488     out->FaceElements->Color[k]=0;
1489    
1490     if( useElementsOnFace )
1491     {
1492     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1493     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1494     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;
1495     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1496     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1497     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1498     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1499     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1500     } else {
1501     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1502     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1503     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1504     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1505     }
1506     }
1507     faceNECount += NE2;
1508     totalNECount += NE2;
1509     }
1510    
1511     }
1512     if( !periodic[2] ) {
1513     /* x3=0 */
1514     for( i1=0; i1<NE1; i1++ )
1515     {
1516     k = i1+faceNECount;
1517     node0 = numDOFLocal*(i1+1) - 1;
1518     node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;
1519    
1520     out->FaceElements->Tag[k]=200;
1521     out->FaceElements->Id[k]=idCount++;
1522     out->FaceElements->Color[k]=0;
1523    
1524     if( useElementsOnFace )
1525     {
1526     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1527     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1528     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1529     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1530     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1531     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1532     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1533     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;
1534     } else {
1535     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1536     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1537     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1538     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1539     }
1540     }
1541     faceNECount += NE1;
1542     totalNECount += NE1;
1543    
1544     /* x3=1 */
1545     for( i1=0; i1<NE1; i1++ )
1546     {
1547     k = i1+faceNECount;
1548     node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;
1549     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;
1550    
1551     out->FaceElements->Tag[k]=200;
1552     out->FaceElements->Id[k]=idCount++;
1553     out->FaceElements->Color[k]=0;
1554    
1555     if( useElementsOnFace )
1556     {
1557     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1558     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1559     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;
1560     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;
1561     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1562     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1563     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1564     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1565     } else {
1566     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1567     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1568     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1569     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1570     }
1571     }
1572     faceNECount += NE1;
1573     totalNECount += NE1;
1574     }
1575     if( !periodic[1] ) {
1576     /* x2=0 */
1577     for( i2=0; i2<NE2; i2++ )
1578     {
1579     k = i2+faceNECount;
1580     node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;
1581     node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1582    
1583     out->FaceElements->Tag[k]=20;
1584     out->FaceElements->Id[k]=idCount++;
1585     out->FaceElements->Color[k]=0;
1586    
1587     if( useElementsOnFace )
1588     {
1589     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1590     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1591     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1592     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1593     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1594     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;
1595     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1596     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1597     } else {
1598     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1599     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1600     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1601     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1602     }
1603     }
1604     faceNECount += NE2;
1605     totalNECount += NE2;
1606    
1607     /* x2=1 */
1608     for( i2=0; i2<NE2; i2++ )
1609     {
1610     k = i2+faceNECount;
1611     node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;
1612     node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1613    
1614     out->FaceElements->Tag[k]=20;
1615     out->FaceElements->Id[k]=idCount++;
1616     out->FaceElements->Color[k]=0;
1617    
1618     if( useElementsOnFace ){
1619     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1620     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;
1621     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;
1622     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;
1623     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1624     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1625     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1626     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1627     } else {
1628     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1629     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1630     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1631     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1632     }
1633     }
1634     faceNECount += NE2;
1635     totalNECount += NE2;
1636     }
1637     }
1638     out->FaceElements->minColor=0;
1639     out->FaceElements->maxColor=0;//23;
1640    
1641     out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;
1642     out->FaceElements->elementDistribution->numLocal = faceNECount;
1643    
1644    
1645     /* setup distribution info for other elements */
1646     out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
1647     out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
1648    
1649     /* condense the nodes: */
1650     Finley_Mesh_resolveNodeIds( out );
1651    
1652     /* setup the CommBuffer */
1653     Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1654     if ( !Finley_MPI_noError( mpi_info )) {
1655     if( Finley_noError() )
1656     Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1657     Paso_MPIInfo_dealloc( mpi_info );
1658     Finley_Mesh_dealloc(out);
1659     return NULL;
1660     }
1661    
1662     Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1663    
1664     /* prepare mesh for further calculatuions:*/
1665     Finley_Mesh_prepare(out) ;
1666    
1667     // print_mesh_statistics( out );
1668    
1669     #ifdef Finley_TRACE
1670     printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1671     #endif
1672    
1673     if (! Finley_noError()) {
1674     Finley_Mesh_dealloc(out);
1675     return NULL;
1676     }
1677    
1678     return out;
1679     }
1680     #endif
1681    

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26