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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 759 by bcumming, Thu Jun 29 01:53:23 2006 UTC revision 1811 by ksteube, Thu Sep 25 23:11:13 2008 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2008 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open Software License version 3.0    *  *
8   *     http://www.opensource.org/licenses/osl-3.0.php       *  * Primary Business: Queensland, Australia
9   *                                                          *  * Licensed under the Open Software License version 3.0
10   ************************************************************  * http://www.opensource.org/licenses/osl-3.0.php
11  */  *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
# Line 20  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$ */  
   
 /**************************************************************/  
   
25  #include "RectangularMesh.h"  #include "RectangularMesh.h"
26    
27  /**************************************************************/  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,
28                                              double* Length,
29  #ifdef PASO_MPI                                            bool_t* periodic,
30  /* get the number of nodes/elements for domain with rank=rank, of size processors                                            index_t order,
31     where n is the total number of nodes/elements in the global domain */                                            index_t reduced_order,
32  static index_t domain_MODdim( index_t rank, index_t size, index_t n )                                            bool_t useElementsOnFace,
33  {                                            bool_t useFullElementOrder,
34    rank = size-rank-1;                                            bool_t optimize)
   
   if( rank < n%size )  
     return (index_t)floor(n/size)+1;  
   return (index_t)floor(n/size);  
 }  
   
   
 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */  
 /* A bit messy, but it only has to be done once... */  
 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 )  
35  {  {
36    index_t i0;    #define N_PER_E 1
37    dim_t numNodesGlobal = numElementsGlobal+1;    #define DIM 3
38      dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2, local_N0, local_N1, local_N2;
39    (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );    dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, NN;
40        index_t node0, myRank, e_offset2, e_offset1, e_offset0, offset1, offset2, offset0, global_i0, global_i1, global_i2;
41    numElementsLocal[0] = numNodesLocal[0]+1;    Finley_Mesh* out;
42    periodicLocal[0] = periodicLocal[1] = FALSE;    Paso_MPIInfo *mpi_info = NULL;
43    nodesExternal[0] = nodesExternal[1] = 1;    char name[50];
44    if( periodic )    double time0=Finley_timer();
   {  
     if( size==1 )  
     {  
       numElementsLocal[0] = numElementsGlobal;  
       nodesExternal[0] = nodesExternal[1] = 0;  
       periodicLocal[0] = periodicLocal[1] = TRUE;  
     }  
     else  
     {  
       if( rank==0 )  
       {  
         periodicLocal[0] = TRUE;  
         numNodesLocal[0]++;  
       }  
       if( rank==(size-1) )  
       {        
         periodicLocal[1] = TRUE;  
         numNodesLocal[0]--;  
         numElementsLocal[0]--;  
       }  
     }  
   }  
   else if( !periodic )  
   {  
     if( rank==0 ){  
       nodesExternal[0]--;  
       numElementsLocal[0]--;  
     }  
     if( rank==(size-1) )  
     {  
       nodesExternal[1]--;  
       numElementsLocal[0]--;  
     }  
   }  
   numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];  
     
   numElementsInternal[0] = numElementsLocal[0];  
   if( (rank==0) && (rank==size-1) );  
       
   else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )  
       numElementsInternal[0] -= 1;  
   else  
     numElementsInternal[0] -= 2;  
   
   firstNode[0] = 0;  
   for( i0=0; i0<rank; i0++ )  
     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );  
   
   numDOFLocal[0] = numNodesLocal[0];  
   if( periodicLocal[0] )  
   {  
     numDOFLocal[0]--;  
   }  
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
 }  
   
 void print_mesh_statistics( Finley_Mesh *out )  
 {  
   index_t i, j, N;  
     
   printf( "\nNodes\n=====\n\n" );  
     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);  
   for( i=0; i<out->Nodes->numNodes; i++ )  
     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)] );  
   
   printf( "Elements\n========\n\n" );  
     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 );  
     N = out->Elements->ReferenceElement->Type->numNodes;  
   for( i=0; i<out->Elements->numElements; i++ ){  
     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)] );  
         printf( " DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(0,i,N)]]] );    
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(j,i,N)]]]  );    
         printf( " ]\n" );    
   }  
45    
46      printf( "\nFace Elements\n==============\n\n" );    /* get MPI information */
47      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 );    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
48      N = out->FaceElements->ReferenceElement->Type->numNodes;    if (! Finley_noError()) {
49    for( i=0; i<out->FaceElements->numElements; i++ ){          return NULL;
         printf( "face element %d \t: id %d  \t: nodes [ %3d", i, out->FaceElements->Id[i], out->FaceElements->Nodes[INDEX2(0,i,N)] );  
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->FaceElements->Nodes[INDEX2(j,i,N)]  );      
         printf( " ] DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(0,i,N)]]] );  
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,N)]]]  );    
         printf( " ]\n" );    
50    }    }
51  }    myRank=mpi_info->rank;
52    
53  #endif    /* set up the global dimensions of the mesh */
54    
 #ifndef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  
 #else  
 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  
 #endif  
 {  
   dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;  
   index_t node0;  
   Finley_Mesh* out;  
   char name[50];  
   double time0=Finley_timer();  
55    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
56    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
57    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
58    N0=NE0+1;    N0=N_PER_E*NE0+1;
59    N1=NE1+1;    N1=N_PER_E*NE1+1;
60    N2=NE2+1;    N2=N_PER_E*NE2+1;
 #ifdef PASO_MPI  
   Paso_MPIInfo *mpi_info = NULL;  
61    
62    /* get MPI information */    /*  allocate mesh: */  
63    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
64    if (! Finley_noError())    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65          return NULL;    if (! Finley_noError()) {
66  #endif        Paso_MPIInfo_free( mpi_info );
67          return NULL;
68      }
69    
70    if (N0<=MIN(N1,N2)) {   Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info));
71       if (N1 <= N2) {   if (useElementsOnFace) {
72          M0=1;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex8Face,
73          M1=N0;                                                                    out->order,
74          M2=N0*N1;                                                                    out->reduced_order,
75       } else {                                                                    mpi_info));
76          M0=1;           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex8Face_Contact,
77          M2=N0;                                                                      out->order,
78          M1=N0*N2;                                                                      out->reduced_order,
79       }                                                                      mpi_info));
   } else if (N1<=MIN(N2,N0)) {  
      if (N2 <= N0) {  
         M1=1;  
         M2=N1;  
         M0=N2*N1;  
      } else {  
         M1=1;  
         M0=N1;  
         M2=N1*N0;  
      }  
80    } else {    } else {
81       if (N0 <= N1) {           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4,
82          M2=1;                                                                    out->order,
83          M0=N2;                                                                    out->reduced_order,
84          M1=N2*N0;                                                                    mpi_info));
85       } else {           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4_Contact,
86          M2=1;                                                                       out->order,
87          M1=N2;                                                                       out->reduced_order,
88          M0=N1*N2;                                                                       mpi_info));
89       }    }
90      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
91                                                     out->order,
92                                                     out->reduced_order,
93                                                     mpi_info));
94      if (! Finley_noError()) {
95          Paso_MPIInfo_free( mpi_info );
96          Finley_Mesh_free(out);
97          return NULL;
98    }    }
99    
100      /* work out the largest dimension */
101      if (N2==MAX3(N0,N1,N2)) {
102         Nstride0=1;
103         Nstride1=N0;
104         Nstride2=N0*N1;
105         local_NE0=NE0;
106         e_offset0=0;
107         local_NE1=NE1;
108         e_offset1=0;
109         Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
110      } else if (N1==MAX3(N0,N1,N2)) {
111         Nstride0=N2;
112         Nstride1=N0*N2;
113         Nstride2=1;
114         local_NE0=NE0;
115         e_offset0=0;
116         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
117         local_NE2=NE2;
118         e_offset2=0;
119      } else {
120         Nstride0=N1*N2;
121         Nstride1=1;
122         Nstride2=N1;
123         Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
124         local_NE1=NE1;
125         e_offset1=0;
126         local_NE2=NE2;
127         e_offset2=0;
128      }
129      offset0=e_offset0*N_PER_E;
130      offset1=e_offset1*N_PER_E;
131      offset2=e_offset2*N_PER_E;
132      local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
133      local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
134      local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
135    
136      /* get the number of surface elements */
137    
138    NFaceElements=0;    NFaceElements=0;
139    if (!periodic[0]) {    if (!periodic[2] && (local_NE2>0)) {
140        NDOF0=N0;      NDOF2=N2;
141        NFaceElements+=2*NE1*NE2;      if (offset2==0) NFaceElements+=local_NE1*local_NE0;
142    } else {      if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
       NDOF0=N0-1;  
   }  
   if (!periodic[1]) {  
       NDOF1=N1;  
       NFaceElements+=2*NE0*NE2;  
   } else {  
       NDOF1=N1-1;  
   }  
   if (!periodic[2]) {  
       NDOF2=N2;  
       NFaceElements+=2*NE0*NE1;  
143    } else {    } else {
144        NDOF2=N2-1;        NDOF2=N2-1;
145    }    }
146      
147    /*  allocate mesh: */    if (!periodic[0] && (local_NE0>0)) {
148           NDOF0=N0;
149    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);       if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
150         if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
 #ifndef PASO_MPI  
   out=Finley_Mesh_alloc(name,3,order);  
 #else  
   out=Finley_Mesh_alloc(name,3,order,mpi_info);  
 #endif  
   if (! Finley_noError()) return NULL;  
   
 #ifdef PASO_MPI  
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);  
151    } else {    } else {
152       out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);        NDOF0=N0-1;
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);  
153    }    }
154    out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);    if (!periodic[1] && (local_NE1>0)) {
155  #else       NDOF1=N1;
156    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);       if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
157    if (useElementsOnFace) {       if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);  
158    } else {    } else {
159       out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);        NDOF1=N1-1;
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order);  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);  
 #ifdef PASO_MPI  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );  
 #endif  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
160    }    }
     
   /*  set nodes: */  
161    
162    #pragma omp parallel for private(i0,i1,i2,k)    /*  allocate tables: */
   for (i2=0;i2<N2;i2++) {  
     for (i1=0;i1<N1;i1++) {  
       for (i0=0;i0<N0;i0++) {  
         k=M0*i0+M1*i1+M2*i2;  
         out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];  
         out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
         out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
         out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;  
         out->Nodes->Tag[k]=0;  
         out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);  
       }  
     }  
   }  
   /* tags for the faces: */  
   /* tags for the faces: */  
   if (!periodic[2]) {  
      for (i1=0;i1<N1;i1++) {  
        for (i0=0;i0<N0;i0++) {  
          out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;  
          out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;  
        }  
      }  
   }  
   if (!periodic[1]) {  
     for (i2=0;i2<N2;i2++) {  
       for (i0=0;i0<N0;i0++) {  
          out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;  
          out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;  
       }  
     }  
   }  
   if (!periodic[0]) {  
     for (i2=0;i2<N2;i2++) {  
       for (i1=0;i1<N1;i1++) {  
         out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;  
         out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;  
       }  
     }  
   }  
163    
164    /*   set the elements: */    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
165        Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
166    #pragma omp parallel for private(i0,i1,i2,k,node0)    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
   for (i2=0;i2<NE2;i2++) {  
     for (i1=0;i1<NE1;i1++) {  
       for (i0=0;i0<NE0;i0++) {  
         k=i0+NE0*i1+NE0*NE1*i2;  
         node0=i0+i1*N0+N0*N1*i2;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;  
167    
168        }    if (Finley_noError()) {
169      }       /* create nodes */
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);  
     
   /*   face elements: */  
     
   if  (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=4;  
   }  
   totalNECount=NE0*NE1*NE2;  
   faceNECount=0;  
     
   /*   these are the quadrilateral elements on boundary 1 (x3=0): */  
   if (!periodic[2]) {  
      /* **  elements on boundary 100 (x3=0): */  
     
      #pragma omp parallel for private(i0,i1,k,node0)  
      for (i1=0;i1<NE1;i1++) {  
        for (i0=0;i0<NE0;i0++) {  
          k=i0+NE0*i1+faceNECount;  
          node0=i0+i1*N0;  
170        
171           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;       #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
172           out->FaceElements->Tag[k]=100;       for (i2=0;i2<local_N2;i2++) {
173           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);         for (i1=0;i1<local_N1;i1++) {
174             for (i0=0;i0<local_N0;i0++) {
175           if  (useElementsOnFace) {             k=i0+local_N0*i1+local_N0*local_N1*i2;
176              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             global_i0=i0+offset0;
177              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;             global_i1=i1+offset1;
178              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;             global_i2=i2+offset2;
179              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
180              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
181              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;             out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
182              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
183              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;             out->Nodes->Tag[k]=0;
184           } else {             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
185              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                                                 +Nstride1*(global_i1%NDOF1)
186              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;                                                 +Nstride2*(global_i2%NDOF2);
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
187           }           }
188         }         }
189       }       }
190       totalNECount+=NE1*NE0;       /*   set the elements: */
191       faceNECount+=NE1*NE0;       NN=out->Elements->numNodes;
192           #pragma omp parallel for private(i0,i1,i2,k,node0)
193       /* **  elements on boundary 200 (x3=1) */       for (i2=0;i2<local_NE2;i2++) {
194             for (i1=0;i1<local_NE1;i1++) {
195       #pragma omp parallel for private(i0,i1,k,node0)           for (i0=0;i0<local_NE0;i0++) {
196       for (i1=0;i1<NE1;i1++) {            
197         for (i0=0;i0<NE0;i0++) {             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
198           k=i0+NE0*i1+faceNECount;             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
199           node0=i0+i1*N0+N0*N1*(NE2-1);    
200                 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
201           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;             out->Elements->Tag[k]=0;
202           out->FaceElements->Tag[k]=200;             out->Elements->Owner[k]=myRank;
203           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;  
204               out->Elements->Nodes[INDEX2(0,k,NN)] =node0                             ;
205           if  (useElementsOnFace) {             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                    Nstride0;
206              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;             out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
207              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;             out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1;
208              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                   ;
209              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;             out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2         +Nstride0;
210              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
211              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1           ;
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;  
212           }           }
   
     
213         }         }
214       }       }
215       totalNECount+=NE1*NE0;       /* face elements */
216       faceNECount+=NE1*NE0;       NN=out->FaceElements->numNodes;
217    }       totalNECount=NE0*NE1*NE2;
218    if (!periodic[0]) {       faceNECount=0;
219       /* **  elements on boundary 001 (x1=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
220           if (!periodic[2]  && (local_NE2>0)) {
221       #pragma omp parallel for private(i1,i2,k,node0)         /* **  elements on boundary 100 (x3=0): */
222       for (i2=0;i2<NE2;i2++) {         if (e_offset2==0) {
223         for (i1=0;i1<NE1;i1++) {            #pragma omp parallel for private(i0,i1,k,node0)
224           k=i1+NE1*i2+faceNECount;            for (i1=0;i1<local_NE1;i1++) {
225           node0=i1*N0+N0*N1*i2;              for (i0=0;i0<local_NE0;i0++) {
226                
227           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;                k=i0+local_NE0*i1+faceNECount;
228           out->FaceElements->Tag[k]=1;                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
          out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;  
          }  
        }  
      }  
      totalNECount+=NE1*NE2;  
      faceNECount+=NE1*NE2;  
229            
230       /* **  elements on boundary 002 (x1=1): */                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
231                    out->FaceElements->Tag[k]=100;
232       #pragma omp parallel for private(i1,i2,k,node0)                out->FaceElements->Owner[k]=myRank;
233       for (i2=0;i2<NE2;i2++) {          
234         for (i1=0;i1<NE1;i1++) {                if  (useElementsOnFace) {
235           k=i1+NE1*i2+faceNECount;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
236           node0=(NE0-1)+i1*N0+N0*N1*i2 ;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0         +Nstride1           ;
237                       out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0         +Nstride1+Nstride0;
238           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+          Nstride0           ;
239           out->FaceElements->Tag[k]=2;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                      ;
240           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1           ;
241                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
242           if  (useElementsOnFace) {                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2          +Nstride0;
243              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                } else {
244              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
245              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           Nstride1           ;
246              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
247              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                     Nstride0;
248              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;                }
249              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;              }
250              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;            }
251           } else {            faceNECount+=local_NE1*local_NE0;
252              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;         }
253              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;         totalNECount+=NE1*NE0;
254              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;         /* **  elements on boundary 200 (x3=1) */
255              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;         if (local_NE2+e_offset2 == NE2) {
256           }            #pragma omp parallel for private(i0,i1,k,node0)
257              for (i1=0;i1<local_NE1;i1++) {
258                for (i0=0;i0<local_NE0;i0++) {
259          
260                  k=i0+local_NE0*i1+faceNECount;
261                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
262            
263                  out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
264                  out->FaceElements->Tag[k]=200;
265                  out->FaceElements->Owner[k]=myRank;
266                  if  (useElementsOnFace) {
267                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
268                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+           Nstride0;
269                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
270                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
271          
272                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
273                     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0                      ;
274                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           Nstride1+Nstride0;
275                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           Nstride1;
276                  } else {
277                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
278                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2           +Nstride0;
279                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
280                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
281                  }
282                }
283              }
284              faceNECount+=local_NE1*local_NE0;
285         }         }
286           totalNECount+=NE1*NE0;
287       }       }
288       totalNECount+=NE1*NE2;       if (!periodic[0]  && (local_NE0>0)) {
289       faceNECount+=NE1*NE2;          /* **  elements on boundary 001 (x1=0): */
290    }      
291    if (!periodic[1]) {          if (e_offset0 == 0) {
292       /* **  elements on boundary 010 (x2=0): */             #pragma omp parallel for private(i1,i2,k,node0)
293                 for (i2=0;i2<local_NE2;i2++) {
294       #pragma omp parallel for private(i0,i2,k,node0)               for (i1=0;i1<local_NE1;i1++) {
295       for (i2=0;i2<NE2;i2++) {        
296         for (i0=0;i0<NE0;i0++) {                 k=i1+local_NE1*i2+faceNECount;
297           k=i0+NE0*i2+faceNECount;                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
298           node0=i0+N0*N1*i2;                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
299                     out->FaceElements->Tag[k]=1;
300           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 out->FaceElements->Owner[k]=myRank;
301           out->FaceElements->Tag[k]=10;        
302           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;                 if  (useElementsOnFace) {
303                      out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
304           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
305              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
306              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1                      ;
307              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0                      ;
308              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;                    out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0           ;
309              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
310              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0           ;
311              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;                 } else {
312              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
313           } else {                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
314              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
315              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1           ;
316              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;                 }
317              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;               }
318           }             }
319         }             faceNECount+=local_NE1*local_NE2;
320       }          }
321       totalNECount+=NE0*NE2;          totalNECount+=NE1*NE2;
322       faceNECount+=NE0*NE2;          /* **  elements on boundary 002 (x1=1): */
323              if (local_NE0+e_offset0 == NE0) {
324       /* **  elements on boundary 020 (x2=1): */             #pragma omp parallel for private(i1,i2,k,node0)
325                 for (i2=0;i2<local_NE2;i2++) {
326       #pragma omp parallel for private(i0,i2,k,node0)               for (i1=0;i1<local_NE1;i1++) {
327       for (i2=0;i2<NE2;i2++) {                 k=i1+local_NE1*i2+faceNECount;
328         for (i0=0;i0<NE0;i0++) {                 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
329           k=i0+NE0*i2+faceNECount;                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
330           node0=i0+(NE1-1)*N0+N0*N1*i2;                 out->FaceElements->Tag[k]=2;
331                     out->FaceElements->Owner[k]=myRank;
332           out->FaceElements->Tag[k]=20;    
333           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 if  (useElementsOnFace) {
334           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      Nstride0;
335                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
336           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
337              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
338              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;        
339              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
340              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           Nstride1           ;
341              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1           ;
342              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2                      ;
343              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;        
344              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;                 } else {
345           } else {                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +Nstride0;
346              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
347              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
348              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
349              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;                 }
350            
351                 }
352               }
353               faceNECount+=local_NE1*local_NE2;
354           }           }
355             totalNECount+=NE1*NE2;
        }  
      }  
      totalNECount+=NE0*NE2;  
      faceNECount+=NE0*NE2;  
   }  
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=23;  
     
   /*  face elements done: */  
     
 #ifdef PASO_MPI  
   /* make sure that the trivial distribution data is correct */  
     out->FaceElements->elementDistribution->numBoundary = 0;  
   out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = faceNECount;  
     out->Elements->elementDistribution->numBoundary = 0;  
   out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal = out->Elements->numElements;  
   out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;  
   out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;  
       
     out->Nodes->degreeOfFreedomDistribution->numInternal = out->Nodes->degreeOfFreedomDistribution->numLocal;  
   out->Nodes->degreeOfFreedomDistribution->numBoundary = 0;  
 #endif  
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
   
 #ifdef PASO_MPI  
   /* setup the CommBuffer */  
   Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
   if ( !Finley_MPI_noError( mpi_info )) {  
     if( Finley_noError() )  
       Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
 #endif  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
   
   #ifdef Finley_TRACE  
   printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);  
   #endif  
   
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   return out;  
 }  
   
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  
 {  
   dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;  
   dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;  
   bool_t dom_left, dom_right, dom_internal;  
   
   index_t N0dom;  
   index_t firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1, node2;  
   index_t targetDomain=-1;  
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;  
   index_t *indexForward=NULL;  
   Finley_Mesh* out;  
   
   char name[50];  
   Paso_MPIInfo *mpi_info = NULL;  
   double time0=Finley_timer();  
   
   NE0=MAX(1,numElements[0]);  
   NE1=MAX(1,numElements[1]);  
   NE2=MAX(1,numElements[2]);  
   N0=NE0+1;  
   N1=NE1+1;  
   N2=NE2+1;  
   
   
   /* get MPI information */  
   mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  
   if (! Finley_noError())  
         return NULL;  
   
   /* use the serial version to generate the mesh for the 1-CPU case */  
   if( mpi_info->size==1 )  
   {  
     Paso_MPIInfo_dealloc( mpi_info );  
     out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace );  
     //print_mesh_statistics( out );  
         return out;  
   }      
   
   if( mpi_info->rank==0 )  
     domLeft = TRUE;  
   if( mpi_info->rank==mpi_info->size-1 )  
     domRight = TRUE;  
   if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )  
     domInternal = TRUE;  
   
   /* dimensions of the local subdomain */  
   domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );    
   
   /* count Degrees of Freedom along each axis */  
   NDOF0 = N0 - periodic[0];  
   NDOF1 = N1 - periodic[1];  
   NDOF2 = N2 - periodic[2];  
   
   /* count face elements */  
   /* internal face elements */  
   NFaceElements = 0;  
   if( !periodic[0] )  
     NFaceElements += (domLeft+domRight)*NE1*NE2;  
   if( !periodic[1] )  
     NFaceElements += 2*numElementsInternal*NE2;  
   if( !periodic[2] )  
     NFaceElements += 2*numElementsInternal*NE1;  
   /* boundary face elements */  
     /* this is looks nasty, but it beats a bunch of nested if/then/else carry-on */  
   NFaceElements += 2*( 2 - (domLeft + domRight)*(!periodic[0]) )*( (!periodic[1])*NE2 + (!periodic[2])*NE1 );  
   
     
   /*  allocate mesh: */  
   sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);  
   
   out=Finley_Mesh_alloc(name,3,order,mpi_info);  
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);  
   
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,(numNodesLocal+2-!periodic[0]*(domLeft+domRight))*N1*N2);  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, (DOFExternal[0]+DOFExternal[1])*NDOF1*NDOF2, 0 );  
   Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
     
   /*  set nodes: */  
   /* INTERNAL/BOUNDARY NODES */  
     k=0;  
   #pragma omp parallel for private(i0,i1,i2,k)  
   for (i2=0;i2<N2;i2++) {  
     for (i1=0;i1<N1;i1++) {  
       for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {          
         out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];  
         out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
         out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
         out->Nodes->Id[k]=k;  
         out->Nodes->Tag[k]=0;  
         out->Nodes->degreeOfFreedom[k]=k;  
       }  
     }  
   }  
   if( domLeft && periodic[0] ) {  
     for (i2=0;i2<N2;i2++)  {  
       for (i1=0;i1<N1;i1++, k++) {  
         out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];  
         out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
         out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
         out->Nodes->Id[k]=k;  
         out->Nodes->Tag[k]=0;  
         out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;  
       }  
     }  
     /* tag the faces for this special case */  
     if( !periodic[1] )  
     {  
       for( i2=0; i2<N2; i2++ ){  
         out->Nodes->Tag[k + (i2-N2)*N1     ] += 10;  
         out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;  
       }  
     }  
     if( !periodic[2] )  
     {  
       for( i1=0; i1<N1; i1++ ){  
         out->Nodes->Tag[k -N1*N2 +i1] += 100;  
         out->Nodes->Tag[k -N1    +i1] += 200;        
       }  
     }  
   }  
   /* tags for the faces: */  
   N0dom = (numNodesLocal-periodicLocal[0]);  
   if (!periodic[2]) {      
     for (i1=0;i1<N1;i1++) {  
       for (i0=0;i0<N0dom;i0++) {    
          out->Nodes->Tag[i0 + N0dom*i1]+=100;  
          out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;  
        }  
356       }       }
357    }       if (!periodic[1] && (local_NE1>0)) {
358    if (!periodic[1]) {          /* **  elements on boundary 010 (x2=0): */
359      for (i2=0;i2<N2;i2++) {          if (e_offset1 == 0) {
360        for (i0=0;i0<N0dom;i0++) {             #pragma omp parallel for private(i0,i2,k,node0)
361           out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;             for (i2=0;i2<local_NE2;i2++) {
362           out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;               for (i0=0;i0<local_NE0;i0++) {
363        }                 k=i0+local_NE0*i2+faceNECount;
364      }                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
365    }          
366    if (!periodic[0] && !domInternal ) {                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
367      for (i2=0;i2<N2;i2++) {                 out->FaceElements->Tag[k]=10;
368        for (i1=0;i1<N1;i1++) {                 out->FaceElements->Owner[k]=myRank;
369          if( domLeft )                 if  (useElementsOnFace) {
370            out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
371          if( domRight )                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
372            out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2           +Nstride0;
373        }                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
374      }        
375    }                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           Nstride1           ;
376      /* setup the forward communication data for the boundary nodes that we have just defined */                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+           Nstride0;
377      /* the case where there are 2 subdomains and periodic[0]=true has to be treated                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
378           as a special case to because the two domains have two interface boundaries to one-another */                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1           ;
379    indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );                 } else {
380      if( mpi_info->size>2 || !periodic[0] ){                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
381          if( domInternal || domRight || periodicLocal[0] )                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
382          {                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+           Nstride0;
383                  for( int i=0; i<NDOF2; i++ )                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
384                      for( int j=0; j<NDOF1; j++ )                 }
385                          indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;               }
386                  targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;             }
387                  Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );             faceNECount+=local_NE0*local_NE2;
         }  
         if( domInternal || domLeft || periodicLocal[1] )  
         {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;  
             targetDomain = (mpi_info->rank+1) % mpi_info->size;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
         }  
     }  
     else {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;  
             targetDomain = (mpi_info->rank+1) % mpi_info->size;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
   
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;  
             targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
     }  
       
   /* EXTERNAL NODES */  
   /* left hand boundary */  
   DOFcount = NDOF1*NDOF2*numDOFLocal;  
   if( (domLeft && periodic[0]) || !domLeft ) {  
     if( (domLeft && periodic[0]) )  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }  
     else  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
388          }          }
389        }                totalNECount+=NE0*NE2;
390          DOFcount += NDOF1*NDOF2;          /* **  elements on boundary 020 (x2=1): */
391          targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;          if (local_NE1+e_offset1 == NE1) {
392          if( !periodic[1] ){             #pragma omp parallel for private(i0,i2,k,node0)
393              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );             for (i2=0;i2<local_NE2;i2++) {
394          }               for (i0=0;i0<local_NE0;i0++) {
395          else {                 k=i0+local_NE0*i2+faceNECount;
396              for( int i=0; i<NDOF2; i++ )                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
397                  for( int j=0; j<NDOF1; j++ )    
398                      indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
399              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );                 out->FaceElements->Tag[k]=20;
400          }                 out->FaceElements->Owner[k]=myRank;
401                          
402      /* tag the faces for this special case */                 if  (useElementsOnFace) {
403      if( !periodic[1] )                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
404      {                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
405        for( i1=0; i1<N1; i1++ ){                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
406          out->Nodes->Tag[k -N1*N2 +i1] += 10;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0           ;
407          out->Nodes->Tag[k -N1    +i1] += 20;        
408        }                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
409      }                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2                      ;
410      if( periodic[2] )                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+           Nstride0;
411      {                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      Nstride0;
412        for( i2=0; i2<N2; i2++ ){                 } else {
413          out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
414          out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;                          out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
415        }                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
416      }                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           Nstride1+Nstride0;
417    }                 }
418    if( (domRight && periodic[0]) || !domRight )               }
419    {             }
420      if( domRight && periodic[0] )             faceNECount+=local_NE0*local_NE2;
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
421          }          }
422        }          totalNECount+=NE0*NE2;
     else  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }  
         DOFcount += NDOF1*NDOF2;  
   
         targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;  
         if( !periodic[1] ){  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );  
         }  
         else {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
         }  
               
     /* tag the faces for this special case */  
     if( !periodic[1] )  
     {  
       for( i1=0; i1<N1; i1++ ){  
         out->Nodes->Tag[k -N1*N2 +i1] += 10;  
         out->Nodes->Tag[k -N1    +i1] += 20;  
       }  
     }  
     if( !periodic[2] )  
     {  
       for( i2=0; i2<N2; i2++ ){  
         out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;  
         out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;        
       }  
     }  
   }  
     out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));  
     out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;  
     
     TMPMEMFREE( indexForward );  
   /*   set the elements: */  
   
   /* INTERNAL elements */  
   N0dom = (numNodesLocal-periodicLocal[0]);  
   k = 0;  
   #pragma omp parallel for private(i0,i1,i2,k,node0)  
   for (i2=0;i2<NE2;i2++) {  
     for (i1=0;i1<NE1;i1++) {  
       for (i0=0;i0<numElementsInternal;i0++,k++) {  
         node0=i0+i1*N0dom+N0dom*N1*i2;  
   
         out->Elements->Id[k]=k;  
                   
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;  
   
       }  
     }  
   }  
     out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;  
     out->Elements->elementDistribution->numBoundary = 0;  
   
   /* BOUNDARY Elements */  
   /* left boundary */  
   if( !domLeft )  
   {  
     for (i2=0;i2<NE2;i2++) {  
       node0 = numNodesLocal*N1*N2 + i2*N1;  
       for (i1=0;i1<NE1;i1++,node0++,k++) {  
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;      
       }  
     }  
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
   }  
   /* the left periodic boundary is done a little differently to a left internal boundary */  
   else if( (domLeft && periodic[0]) )  
   {  
     for (i2=0;i2<NE2;i2++) {  
       node0 = numDOFLocal*N1*N2 + i2*N1;  
       node1 = node0 + N1*N2;  
       for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;      
       }  
     }    
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
   }  
   /* right boundary */  
   if( !domRight || (domRight && periodic[0]) ){  
     for (i2=0;i2<NE2;i2++) {  
       for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {  
                 node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;  
                 node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;    
       }  
     }  
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
     }  
   
     out->Elements->minColor=0;  
   out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);  
     out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;  
       
   /*   face elements: */  
     
   if  (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=4;  
   }  
   totalNECount=k;  
   faceNECount=0;  
     idCount = totalNECount;  
     
     /*   Do internal face elements for each boundary face first */  
   /*   these are the quadrilateral elements on boundary 1 (x3=0): */  
     numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];  
   if (!periodic[2]) {  
      /*  elements on boundary 100 (x3=0): */  
     
      #pragma omp parallel for private(i0,i1,k,node0)  
      for (i1=0;i1<NE1;i1++) {  
        for (i0=0; i0<numElementsInternal; i0++) {  
          k=i0+numElementsInternal*i1+faceNECount;  
          node0=i0+i1*numDOFLocal;  
     
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Tag[k]=100;  
          out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
          }  
        }  
423       }       }
424       totalNECount+=NE1*numElementsInternal;       /* add tag names */
425       faceNECount+=NE1*numElementsInternal;       Finley_Mesh_addTagMap(out,"top", 200);
426                 Finley_Mesh_addTagMap(out,"bottom", 100);
427       /* **  elements on boundary 200 (x3=1) */       Finley_Mesh_addTagMap(out,"left", 1);
428           Finley_Mesh_addTagMap(out,"right", 2);
429       #pragma omp parallel for private(i0,i1,k,node0)       Finley_Mesh_addTagMap(out,"front", 10);
430       for (i1=0;i1<NE1;i1++) {       Finley_Mesh_addTagMap(out,"back", 20);
        for (i0=0;i0<numElementsInternal;i0++) {  
          k=i0+numElementsInternal*i1+faceNECount;  
          node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);  
431        
432           out->FaceElements->Id[k]=idCount++;       /* prepare mesh for further calculatuions:*/
433           out->FaceElements->Tag[k]=200;       if (Finley_noError()) {
434           out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;           Finley_Mesh_resolveNodeIds(out);
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;  
          }  
        }  
      }  
      totalNECount+=NE1*numElementsInternal;  
      faceNECount+=NE1*numElementsInternal;  
   }  
   if (!periodic[0] && !domInternal) {  
      /* **  elements on boundary 001 (x1=0): */  
      if( domLeft ){  
              #pragma omp parallel for private(i1,i2,k,node0)  
              for (i2=0;i2<NE2;i2++) {  
                  for (i1=0;i1<NE1;i1++) {  
                      k=i1+NE1*i2+faceNECount;  
                      node0=i1*numDOFLocal+numDOFLocal*N1*i2;  
           
                      out->FaceElements->Id[k]=idCount++;  
                      out->FaceElements->Tag[k]=1;  
                      out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;  
   
                      if  (useElementsOnFace) {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
                             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;  
                      } else {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;  
                      }  
                  }  
              }  
              totalNECount+=NE1*NE2;  
              faceNECount+=NE1*NE2;  
435       }       }
436       /* **  elements on boundary 002 (x1=1): */       if (Finley_noError()) {
437           if( domRight ) {           Finley_Mesh_prepare(out, optimize);
              #pragma omp parallel for private(i1,i2,k,node0)  
              for (i2=0;i2<NE2;i2++) {  
                  for (i1=0;i1<NE1;i1++) {  
                      k=i1+NE1*i2+faceNECount;  
                      node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;  
           
                      out->FaceElements->Id[k]=idCount++;  
                      out->FaceElements->Tag[k]=2;  
                      out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;  
   
                      if  (useElementsOnFace) {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;  
                      } else {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                      }  
                  }  
              }  
        totalNECount+=NE1*NE2;  
        faceNECount+=NE1*NE2;  
438       }       }
439    }    }
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x2=0): */  
     
      #pragma omp parallel for private(i0,i2,k,node0)  
      for (i2=0;i2<NE2;i2++) {  
        for (i0=0;i0<numElementsInternal;i0++) {  
          k=i0+numElementsInternal*i2+faceNECount;  
          node0=i0+numDOFLocal*N1*i2;  
     
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Tag[k]=10;  
          out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
          }  
        }  
      }  
      totalNECount+=numElementsInternal*NE2;  
      faceNECount+=numElementsInternal*NE2;  
     
      /* **  elements on boundary 020 (x2=1): */  
     
      #pragma omp parallel for private(i0,i2,k,node0)  
      for (i2=0;i2<NE2;i2++) {  
        for (i0=0;i0<numElementsInternal;i0++) {  
          k=i0+numElementsInternal*i2+faceNECount;  
          node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;  
     
          out->FaceElements->Tag[k]=20;  
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;  
          }  
440    
441         }    if (!Finley_noError()) {
442       }        Finley_Mesh_free(out);
      totalNECount+=numElementsInternal*NE2;  
      faceNECount+=numElementsInternal*NE2;  
   }  
     out->FaceElements->elementDistribution->numInternal = faceNECount;  
       
     /* now do the boundary face elements */  
     /* LHS */  
     if( !domLeft  )  
     {  
         if( !periodic[2] ) {  
   
             /* x3=0 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = i1*numNodesLocal;  
                 node1 = numNodesLocal*N1*N2 + i1;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
   
             /* x3=1 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;  
                 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
         }  
   
         if( !periodic[1] ) {  
             /* x2=0 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = i2*numNodesLocal*N1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
   
             /* x2=1 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
         }  
     }  
     
     /* RHS */  
     if( !domRight || periodicLocal[1] )  
     {  
         /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */  
         if( domLeft && periodic[0] ){  
             if( !periodic[2] ) {  
   
                 /* x3=0 */  
                 for( i1=0; i1<NE1; i1++ )  
                 {  
                     k = i1+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i1;  
                     node1 = numNodesLocal*N1*N2 + i1;  
   
                     out->FaceElements->Tag[k]=200;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                     }  
                 }    
                 faceNECount  += NE1;  
                 totalNECount += NE1;  
   
                 /* x3=1 */  
                 for( i1=0; i1<NE1; i1++ )  
                 {  
                     k = i1+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;  
                     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;  
   
                     out->FaceElements->Tag[k]=200;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
                       
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;  
                     }  
                 }    
                 faceNECount  += NE1;  
                 totalNECount += NE1;  
             }  
   
             if( !periodic[1] ) {  
                 /* x2=0 */  
                 for( i2=0; i2<NE2; i2++ )  
                 {  
                     k = i2+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i2*N1;  
                     node1 = numNodesLocal*N1*N2 + i2*N1;  
   
                     out->FaceElements->Tag[k]=20;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                     }  
                 }    
                 faceNECount  += NE2;  
                 totalNECount += NE2;  
   
                 /* x2=1 */  
                 for( i2=0; i2<NE2; i2++ )  
                 {  
                     k = i2+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;  
                     node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;  
   
                     out->FaceElements->Tag[k]=20;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                     }  
                 }    
                 faceNECount  += NE2;  
                 totalNECount += NE2;  
             }  
               
         }  
         if( !periodic[2] ) {  
             /* x3=0 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numDOFLocal*(i1+1) - 1;  
                 node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
   
             /* x3=1 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
                   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
         }  
         if( !periodic[1] ) {  
             /* x2=0 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
   
             /* x2=1 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
                   
                 if( useElementsOnFace ){  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 }  
             }  
             faceNECount  += NE2;  
             totalNECount += NE2;  
         }  
     }  
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=0;//23;  
   
     out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;  
   out->FaceElements->elementDistribution->numLocal = faceNECount;  
   
   
   /* setup distribution info for other elements */  
   out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;  
   out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;  
       
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds( out );  
   
   /* setup the CommBuffer */  
   Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
   if ( !Finley_MPI_noError( mpi_info )) {  
     if( Finley_noError() )  
       Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
443    }    }
444      /* free up memory */
445    Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );    Paso_MPIInfo_free( mpi_info );  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
   
 //  print_mesh_statistics( out );  
   
446    #ifdef Finley_TRACE    #ifdef Finley_TRACE
447    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
448    #endif    #endif
       
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
449    
450    return out;    return out;
451  }  }
 #endif  
   

Legend:
Removed from v.759  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26