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

Diff of /trunk-mpi-branch/finley/src/Mesh_hex8.c

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

revision 1222 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 27  Line 27 
27    
28  #include "RectangularMesh.h"  #include "RectangularMesh.h"
29    
30  /**************************************************************/  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,
31                                              double* Length,
32  #ifdef PASO_MPI                                            bool_t* periodic,
33  /* get the number of nodes/elements for domain with rank=rank, of size processors                                            index_t order,
34     where n is the total number of nodes/elements in the global domain */                                            index_t reduced_order,
35  static index_t domain_MODdim( index_t rank, index_t size, index_t n )                                            bool_t useElementsOnFace,
36  {                                            bool_t useFullElementOrder)
   rank = size-rank-1;  
   
   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 )  
37  {  {
38    index_t i0;    #define N_PER_E 1
39    dim_t numNodesGlobal = numElementsGlobal+1;    #define DIM 3
40      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;
41      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, NN;
42      index_t node0, myRank, e_offset2, e_offset1, e_offset0, offset1, offset2, offset0, global_i0, global_i1, global_i2;
43      Finley_Mesh* out;
44      Paso_MPIInfo *mpi_info = NULL;
45      char name[50];
46      double time0=Finley_timer();
47    
48    (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );    /* get MPI information */
49        mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50    numElementsLocal[0] = numNodesLocal[0]+1;    if (! Finley_noError()) {
51    periodicLocal[0] = periodicLocal[1] = FALSE;          return NULL;
   nodesExternal[0] = nodesExternal[1] = 1;  
   if( periodic )  
   {  
     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]--;  
     }  
52    }    }
53    numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];    myRank=mpi_info->rank;
     
   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;  
54    
55    firstNode[0] = 0;    /*  allocate mesh: */  
56    for( i0=0; i0<rank; i0++ )    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
57      firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
58      if (! Finley_noError()) {
59          Paso_MPIInfo_free( mpi_info );
60          return NULL;
61      }
62    
63    numDOFLocal[0] = numNodesLocal[0];   Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info));
64    if( periodicLocal[0] )   if (useElementsOnFace) {
65    {           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex8Face,
66      numDOFLocal[0]--;                                                                    out->order,
67                                                                      out->reduced_order,
68                                                                      mpi_info));
69             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex8Face_Contact,
70                                                                        out->order,
71                                                                        out->reduced_order,
72                                                                        mpi_info));
73      } else {
74             Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4,
75                                                                      out->order,
76                                                                      out->reduced_order,
77                                                                      mpi_info));
78             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4_Contact,
79                                                                         out->order,
80                                                                         out->reduced_order,
81                                                                         mpi_info));
82      }
83      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
84                                                     out->order,
85                                                     out->reduced_order,
86                                                     mpi_info));
87      if (! Finley_noError()) {
88          Paso_MPIInfo_free( mpi_info );
89          Finley_Mesh_free(out);
90          return NULL;
91    }    }
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
 }  
92    
93  #endif    /* set up the global dimensions of the mesh */
94    
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)  
 #else  
 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_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();  
95    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
96    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
97    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
98    N0=NE0+1;    N0=N_PER_E*NE0+1;
99    N1=NE1+1;    N1=N_PER_E*NE1+1;
100    N2=NE2+1;    N2=N_PER_E*NE2+1;
101    
102    if (N0<=MIN(N1,N2)) {    /* work out the largest dimension */
103       if (N1 <= N2) {    if (N2=MAX3(N0,N1,N2)) {
104          M0=1;       Nstride0=1;
105          M1=N0;       Nstride1=N0;
106          M2=N0*N1;       Nstride2=N0*N1;
107       } else {       local_NE0=NE0;
108          M0=1;       e_offset0=0;
109          M2=N0;       local_NE1=NE1;
110          M1=N0*N2;       e_offset1=0;
111       }       Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
112    } else if (N1<=MIN(N2,N0)) {    } else if (N1=MAX3(N0,N1,N2)) {
113       if (N2 <= N0) {       Nstride0=N2;
114          M1=1;       Nstride1=N0*N2;
115          M2=N1;       Nstride2=1;
116          M0=N2*N1;       local_NE0=NE0;
117       } else {       e_offset0=0;
118          M1=1;       Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
119          M0=N1;       local_NE2=NE2;
120          M2=N1*N0;       e_offset2=0;
      }  
121    } else {    } else {
122       if (N0 <= N1) {       Nstride0=N1*N2;
123          M2=1;       Nstride1=1;
124          M0=N2;       Nstride2=N1;
125          M1=N2*N0;       Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
126       } else {       local_NE1=NE1;
127          M2=1;       e_offset1=0;
128          M1=N2;       local_NE2=NE2;
129          M0=N1*N2;       e_offset2=0;
130       }    }
131    }    offset0=e_offset0*N_PER_E;
132      offset1=e_offset1*N_PER_E;
133      offset2=e_offset2*N_PER_E;
134      local_N0=local_NE0*N_PER_E+1;
135      local_N1=local_NE1*N_PER_E+1;
136      local_N2=local_NE2*N_PER_E+1;
137    
138      /* get the number of surface elements */
139    
140    NFaceElements=0;    NFaceElements=0;
141      if (!periodic[2]) {
142        NDOF2=N2;
143        if (offset2==0) NFaceElements+=local_NE1*local_NE0;
144        if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
145      } else {
146          NDOF2=N2-1;
147    
148    if (!periodic[0]) {    if (!periodic[0]) {
149        NDOF0=N0;       NDOF0=N0;
150        NFaceElements+=2*NE1*NE2;       if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
151         if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
152    } else {    } else {
153        NDOF0=N0-1;        NDOF0=N0-1;
154    }    }
155    if (!periodic[1]) {    if (!periodic[1]) {
156        NDOF1=N1;       NDOF1=N1;
157        NFaceElements+=2*NE0*NE2;       if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
158         if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
159    } else {    } else {
160        NDOF1=N1-1;        NDOF1=N1-1;
161    }    }
   if (!periodic[2]) {  
       NDOF2=N2;  
       NFaceElements+=2*NE0*NE1;  
   } else {  
       NDOF2=N2-1;  
   }  
     
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);  
   
 #ifndef PASO_MPI  
   out=Finley_Mesh_alloc(name,3,order,reduced_order);  
 #else  
   out=Finley_Mesh_alloc(name,3,order,reduced_order,mpi_info);  
 #endif  
   if (! Finley_noError()) return NULL;  
162    
163  #ifdef PASO_MPI    /*  allocate tables: */
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order,mpi_info);  
 #else  
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order, out->reduced_order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order, out->reduced_order);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order, out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order, out->reduced_order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order);  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
164    
165        Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
166    /*  allocate tables: */    Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
   Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);  
167    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
 #ifdef PASO_MPI  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1*NE2, NE0*NE1*NE2);  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
     
   /*  set nodes: */  
   
   #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<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);  
 #ifdef PASO_MPI  
         out->Nodes->Dom[k]=NODE_INTERNAL;  
 #endif  
       }  
     }  
   }  
   /* 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;  
       }  
     }  
   }  
   
   /*   set the elements: */  
     
   #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<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);;  
 #ifdef PASO_MPI  
         out->Elements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
168    
169          out->Elements->Nodes[INDEX2(0,k,8)]=node0;    if (Finley_noError()) {
170          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;       /* create nodes */
         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;  
   
       }  
     }  
   }  
   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;  
171        
172           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;       #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
173           out->FaceElements->Tag[k]=100;       for (i2=0;i2<local_N2;i2++) {
174           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);         for (i1=0;i1<local_N1;i1++) {
175  #ifdef PASO_MPI           for (i0=0;i0<local_N0;i0++) {
176          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             k=i0+local_N0*i1+local_N0*local_N1*i2;
177  #endif             global_i0=i0+offset0;
178               global_i1=i1+offset1;
179           if  (useElementsOnFace) {             global_i2=i2+offset2;
180              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
181              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
182              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;             out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
183              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
184              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;             out->Nodes->Tag[k]=0;
185              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
186              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;                                                 +Nstride1*(global_i1%NDOF1)
187              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;                                                 +Nstride2*(global_i2%NDOF2);
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
188           }           }
189         }         }
190       }       }
191       totalNECount+=NE1*NE0;       /*   set the elements: */
192       faceNECount+=NE1*NE0;       NN=out->Elements->numNodes;
193           #pragma omp parallel for private(i0,i1,i2,k,node0)
194       /* **  elements on boundary 200 (x3=1) */       for (i2=0;i2<local_NE2;i2++) {
195             for (i1=0;i1<local_NE1;i1++) {
196       #pragma omp parallel for private(i0,i1,k,node0)           for (i0=0;i0<local_NE0;i0++) {
197       for (i1=0;i1<NE1;i1++) {            
198         for (i0=0;i0<NE0;i0++) {             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
199           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);
200           node0=i0+i1*N0+N0*N1*(NE2-1);    
201                 out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
202           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;             out->Elements->Tag[k]=0;
203           out->FaceElements->Tag[k]=200;             out->Elements->Owner[k]=myRank;
204           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;  
205  #ifdef PASO_MPI             out->Elements->Nodes[INDEX2(0,k,NN)] =node0                             ;
206          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                    Nstride0;
207  #endif             out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
208               out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1;
209           if  (useElementsOnFace) {             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                   ;
210              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;             out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2         +Nstride0;
211              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
212              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1           ;
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;  
             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+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;  
213           }           }
   
     
214         }         }
215       }       }
216       totalNECount+=NE1*NE0;       /* face elements */
217       faceNECount+=NE1*NE0;       NN=out->FaceElements->numNodes;
218    }       totalNECount=NE0*NE1*NE2;
219    if (!periodic[0]) {       faceNECount=0;
220       /* **  elements on boundary 001 (x1=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
221           if (!periodic[2]) {
222       #pragma omp parallel for private(i1,i2,k,node0)         /* **  elements on boundary 100 (x3=0): */
223       for (i2=0;i2<NE2;i2++) {         if (offset2==0) {
224         for (i1=0;i1<NE1;i1++) {            #pragma omp parallel for private(i0,i1,k,node0)
225           k=i1+NE1*i2+faceNECount;            for (i1=0;i1<local_NE1;i1++) {
226           node0=i1*N0+N0*N1*i2;              for (i0=0;i0<local_NE0;i0++) {
227                
228           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;                k=i0+local_NE0*i1+faceNECount;
229           out->FaceElements->Tag[k]=1;                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
230           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;      
231  #ifdef PASO_MPI                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
232          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                out->FaceElements->Tag[k]=100;
233  #endif                out->FaceElements->Owner[k]=myRank;
234            
235           if  (useElementsOnFace) {                if  (useElementsOnFace) {
236              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
237              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0         +Nstride1           ;
238              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0         +Nstride1+Nstride0;
239              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+          Nstride0           ;
240              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                      ;
241              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1           ;
242              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
243              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2          +Nstride0;
244           } else {                } else {
245              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
246              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           Nstride1           ;
247              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
248              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                     Nstride0;
249           }                }
250                }
251              }
252              faceNECount+=local_NE1*local_NE0;
253           }
254           totalNECount+=NE1*NE0;
255           /* **  elements on boundary 200 (x3=1) */
256           if (local_NE2+e_offset2 == NE2) {
257              #pragma omp parallel for private(i0,i1,k,node0)
258              for (i1=0;i1<local_NE1;i1++) {
259                for (i0=0;i0<local_NE0;i0++) {
260          
261                  k=i0+local_NE0*i1+faceNECount;
262                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
263            
264                  out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
265                  out->FaceElements->Tag[k]=200;
266                  out->FaceElements->Owner[k]=myRank;
267                  if  (useElementsOnFace) {
268                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
269                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+           Nstride0;
270                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
271                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
272          
273                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
274                     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0                      ;
275                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           Nstride1+Nstride0;
276                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           Nstride1;
277                  } else {
278                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
279                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2           +Nstride0;
280                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
281                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
282                  }
283                }
284              }
285              faceNECount+=local_NE1*local_NE0;
286         }         }
287           totalNECount+=NE1*NE0;
288       }       }
289       totalNECount+=NE1*NE2;       if (!periodic[0]) {
290       faceNECount+=NE1*NE2;          /* **  elements on boundary 001 (x1=0): */
291            
292       /* **  elements on boundary 002 (x1=1): */          if (e_offset0 == 0) {
293                 #pragma omp parallel for private(i1,i2,k,node0)
294       #pragma omp parallel for private(i1,i2,k,node0)             for (i2=0;i2<local_NE2;i2++) {
295       for (i2=0;i2<NE2;i2++) {               for (i1=0;i1<local_NE1;i1++) {
296         for (i1=0;i1<NE1;i1++) {        
297           k=i1+NE1*i2+faceNECount;                 k=i1+local_NE1*i2+faceNECount;
298           node0=(NE0-1)+i1*N0+N0*N1*i2 ;                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
299                     out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
300           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;                 out->FaceElements->Tag[k]=1;
301           out->FaceElements->Tag[k]=2;                 out->FaceElements->Owner[k]=myRank;
302           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;        
303  #ifdef PASO_MPI                 if  (useElementsOnFace) {
304          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
305  #endif                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
306                      out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
307           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1                      ;
308              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0                      ;
309              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0           ;
310              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                    out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
311              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;                    out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0           ;
312              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;                 } else {
313              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
314              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
315              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
316           } else {                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1           ;
317              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                   }
318              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;                 }
319              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;               }
320              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;             }
321               faceNECount+=local_NE1*local_NE2;
322            }
323            totalNECount+=NE1*NE2;
324        
325            /* **  elements on boundary 002 (x1=1): */
326            if (local_NE0+e_offset0 == NE0) {
327               #pragma omp parallel for private(i1,i2,k,node0)
328               for (i2=0;i2<local_NE2;i2++) {
329                 for (i1=0;i1<local_NE1;i1++) {
330                   k=i1+local_NE1*i2+faceNECount;
331                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
332                   out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
333                   out->FaceElements->Tag[k]=2;
334                   out->FaceElements->Owner[k]=myRank;
335      
336                   if  (useElementsOnFace) {
337                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      Nstride0;
338                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
339                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
340                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
341          
342                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
343                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           Nstride1           ;
344                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1           ;
345                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2                      ;
346          
347                   } else {
348                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +Nstride0;
349                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
350                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
351                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
352                   }
353            
354                 }
355               }
356               faceNECount+=local_NE1*local_NE2;
357           }           }
358         }           totalNECount+=NE1*NE2;
359       }       }
360       totalNECount+=NE1*NE2;       if (!periodic[1]) {
361       faceNECount+=NE1*NE2;          /* **  elements on boundary 010 (x2=0): */
362    }          if (e_offset1 == 0) {
363    if (!periodic[1]) {             #pragma omp parallel for private(i0,i2,k,node0)
364       /* **  elements on boundary 010 (x2=0): */             for (i2=0;i2<local_NE2;i2++) {
365                   for (i0=0;i0<local_NE0;i0++) {
366       #pragma omp parallel for private(i0,i2,k,node0)                 k=i0+local_NE0*i2+faceNECount;
367       for (i2=0;i2<NE2;i2++) {                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
368         for (i0=0;i0<NE0;i0++) {          
369           k=i0+NE0*i2+faceNECount;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
370           node0=i0+N0*N1*i2;                 out->FaceElements->Tag[k]=10;
371                     out->FaceElements->Owner[k]=myRank;
372           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 if  (useElementsOnFace) {
373           out->FaceElements->Tag[k]=10;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
374           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
375  #ifdef PASO_MPI                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2           +Nstride0;
376          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
377  #endif        
378                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           Nstride1           ;
379           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+           Nstride0;
380              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
381              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1           ;
382              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;                 } else {
383              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
384              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
385              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+           Nstride0;
386              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
387              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;                 }
388           } else {               }
389              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             }
390              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;             faceNECount+=local_NE0*local_NE2;
391              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;          }
392              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;          totalNECount+=NE0*NE2;
393           }          /* **  elements on boundary 020 (x2=1): */
394         }          if (local_NE1+e_offset1 == NE1) {
395               #pragma omp parallel for private(i0,i2,k,node0)
396               for (i2=0;i2<local_NE2;i2++) {
397                 for (i0=0;i0<local_NE0;i0++) {
398                   k=i0+local_NE0*i2+faceNECount;
399                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
400      
401                   out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
402                   out->FaceElements->Tag[k]=20;
403                   out->FaceElements->Owner[k]=myRank;
404          
405                   if  (useElementsOnFace) {
406                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
407                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
408                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
409                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0           ;
410          
411                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
412                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2                      ;
413                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+           Nstride0;
414                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      Nstride0;
415                   } else {
416                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
417                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
418                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
419                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           Nstride1+Nstride0;
420                   }
421                 }
422               }
423               faceNECount+=local_NE0*local_NE2;
424            }
425            totalNECount+=NE0*NE2;
426         }
427         /* add tag names */
428         Finley_Mesh_addTagMap(out,"top", 200);
429         Finley_Mesh_addTagMap(out,"bottom", 100);
430         Finley_Mesh_addTagMap(out,"left", 1);
431         Finley_Mesh_addTagMap(out,"right", 2);
432         Finley_Mesh_addTagMap(out,"front", 10);
433         Finley_Mesh_addTagMap(out,"back", 20);
434      
435         /* prepare mesh for further calculatuions:*/
436         if (Finley_noError()) {
437             Finley_Mesh_resolveNodeIds(out);
438       }       }
439       totalNECount+=NE0*NE2;       if (Finley_noError()) {
440       faceNECount+=NE0*NE2;           Finley_Mesh_prepare(out);
     
      /* **  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<NE0;i0++) {  
          k=i0+NE0*i2+faceNECount;  
          node0=i0+(NE1-1)*N0+N0*N1*i2;  
     
          out->FaceElements->Tag[k]=20;  
          out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;  
          out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;  
 #ifdef PASO_MPI  
         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;  
          }  
   
        }  
441       }       }
      totalNECount+=NE0*NE2;  
      faceNECount+=NE0*NE2;  
442    }    }
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=23;  
     
 #ifdef PASO_MPI  
     Finley_ElementFile_setDomainFlags( out->Elements );  
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
443    
444      /* reorder the degrees of freedom */    if (!Finley_noError()) {
445      Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );        Finley_Mesh_free(out);
 #endif  
       
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
   
   /* 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()) {  
        if ( ! Finley_Mesh_isPrepared(out)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
       Finley_Mesh_dealloc(out);  
   }  
   return out;  
 }  
   
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace)  
 {  
   dim_t N0,N1,N2,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;  
   dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;  
   bool_t dom_left, dom_right, dom_internal;  
   index_t firstNode=0, DOFcount=0, node0, node1, node2;  
   index_t targetDomain=-1, firstNodeConstruct, j;  
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;  
     index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;  
   Finley_Mesh* out;  
   
   char name[50];  
   Paso_MPIInfo *mpi_info = NULL;  
   double time0=Finley_timer();  
   
     index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};  
     index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};  
     index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};  
     index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};  
     index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};  
     index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};  
   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 )  
   {  
     out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, reduced_order, useElementsOnFace, mpi_info );  
         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*numElementsLocal*NE2;  
   if( !periodic[2] )  
     NFaceElements += 2*numElementsLocal*NE1;  
     
     boundaryLeft = !domLeft || periodicLocal[0];  
     boundaryRight = !domRight || periodicLocal[1];  
     N0t = numNodesLocal + boundaryRight + boundaryLeft;  
     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;  
     firstNodeConstruct = firstNode - boundaryLeft;  
     firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;  
   
   /*  allocate mesh: */  
   sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);  
   
   out=Finley_Mesh_alloc(name,3,order,reduced_order,mpi_info);  
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,out->reduced_order, mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,out->reduced_order, mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,out->reduced_order, mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order, mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,out->reduced_order, mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,out->reduced_order, mpi_info);  
   
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);  
   Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );  
   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<N0t;i0++,k++) {          
         out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/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]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;  
                 out->Nodes->Dom[k]=NODE_INTERNAL;  
       }  
     }  
   }  
   
     /* mark the nodes that reference external and boundary DOF as such */  
     if( boundaryLeft ){  
         for( i1=0; i1<N1; i1++ )  
             for( i2=0; i2<N2; i2++ ) {  
                 out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;  
                 out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;  
             }  
     }  
     if( boundaryRight ){  
         for( i1=0; i1<N1; i1++ )  
             for( i2=0; i2<N2; i2++ ) {  
                 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;  
                 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;  
             }  
     }  
     if( periodicLocal[0] ){  
         for( i1=0; i1<N1; i1++ )  
             for( i2=0; i2<N2; i2++ ) {  
                 out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];  
                 out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;  
             }  
     }  
           
   /* tag Nodes that are referenced by face elements */  
   if (!periodic[2]) {      
     for (i1=0;i1<N1;i1++) {  
       for (i0=0;i0<N0t;i0++) {    
          out->Nodes->Tag[i0 + N0t*i1]+=100;  
          out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;  
        }  
      }  
446    }    }
   if (!periodic[1]) {  
     for (i2=0;i2<N2;i2++) {  
       for (i0=0;i0<N0t;i0++) {  
          out->Nodes->Tag[i0 + i2*N1*N0t]+=10;  
          out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;  
       }  
     }  
   }  
   if (!periodic[0] && !domInternal ) {  
     for (i2=0;i2<N2;i2++) {  
       for (i1=0;i1<N1;i1++) {  
         if( domLeft )  
           out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;  
         if( domRight )  
           out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;  
       }  
     }  
   }  
   
     /* form the boudary communication information */  
     forwardDOF  = MEMALLOC(NDOF1*NDOF2,index_t);  
     backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);  
     if( !(mpi_info->size==2 && periodicLocal[0])){  
         if( boundaryLeft  ) {  
             targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;  
             for( i2=0; i2<NDOF2; i2++ ){  
                 for( i1=0; i1<NDOF1; i1++ ){  
                     forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];  
                     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];  
                 }  
             }  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );  
         }  
         if( boundaryRight ) {  
             targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;  
             for( i2=0; i2<NDOF2; i2++ ){  
                 for( i1=0; i1<NDOF1; i1++ ){  
                     forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];  
                     backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];  
                 }  
             }  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );  
         }  
     } else{  
         /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */  
         targetDomain = 1;  
           
         for( i2=0; i2<NDOF2; i2++ ){  
             for( i1=0; i1<NDOF1; i1++ ){  
                 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];  
                 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];  
             }  
         }  
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );  
   
         for( i2=0; i2<NDOF2; i2++ ){  
             for( i1=0; i1<NDOF1; i1++ ){  
                 forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];  
                 backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];  
             }  
         }  
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );  
     }  
     MEMFREE( forwardDOF );  
     MEMFREE( backwardDOF );  
   /*   set the elements: */  
   
   /* INTERNAL elements */  
   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<numElementsLocal;i0++,k++) {  
                 node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t :  i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];  
                   
                 out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
                 out->Elements->Dom[k]=ELEMENT_INTERNAL;  
   
         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+N0t+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;  
   
       }  
     }  
   }  
     out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);  
     if( boundaryLeft )  
         for( i2=0; i2<NE2; i2++ )  
             for( i1=0; i1<NE1; i1++ )  
                 out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;  
     if( boundaryRight )  
         for( i2=0; i2<NE2; i2++ )  
             for( i1=0; i1<NE1; i1++ )  
                 out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
           
     Finley_ElementFile_setDomainFlags( out->Elements );  
       
   /*   face elements: */  
   if  (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=4;  
   }  
   totalNECount=out->Elements->numElements;  
   faceNECount=0;  
     idCount = totalNECount;  
     
   /*   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)  
      for (i1=0;i1<NE1;i1++) {  
        for (i0=0; i0<numElementsLocal; i0++) {  
          k=i0+numElementsLocal*i1+faceNECount;  
                  kk=i0 + i1*numElementsLocal;  
                  facePerm = face2;  
     
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Tag[k]=100;  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);  
   
                  for( j=0; j<NUMNODES; j++ )  
                     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
        }  
      }  
        if( boundaryLeft ){  
             for( i1=0; i1<NE1; i1++ )  
                 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;  
             if( periodicLocal[0] )  
                 for( i1=0; i1<NE1; i1++ )  
                     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;  
          }  
        if( boundaryRight )  
             for( i1=0; i1<NE1; i1++ )  
                 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
      totalNECount+=NE1*numElementsLocal;  
      faceNECount+=NE1*numElementsLocal;  
           
      /* **  elements on boundary 200 (x3=1) */  
     
      #pragma omp parallel for private(i0,i1,k)  
      for (i1=0;i1<NE1;i1++) {  
        for (i0=0;i0<numElementsLocal;i0++) {  
          k=i0+numElementsLocal*i1+faceNECount;  
                  kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);  
                  facePerm = face3;  
     
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Tag[k]=200;  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
          out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;  
   
                  for( j=0; j<NUMNODES; j++ )  
                     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
        }  
      }  
        if( boundaryLeft ){  
             for( i1=0; i1<NE1; i1++ )  
                 out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;  
             if( periodicLocal[0] )  
                 for( i1=0; i1<NE1; i1++ )  
                     out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;  
          }  
        if( boundaryRight )  
             for( i1=0; i1<NE1; i1++ )  
                 out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
      totalNECount+=NE1*numElementsLocal;  
      faceNECount+=NE1*numElementsLocal;  
   }  
   if (!periodic[0] && !domInternal) {  
      /* **  elements on boundary 001 (x1=0): */  
      if( domLeft ){  
              #pragma omp parallel for private(i1,i2,k)  
              for (i2=0;i2<NE2;i2++) {  
                  for (i1=0;i1<NE1;i1++) {  
                      k=i1+NE1*i2+faceNECount;  
                      kk=i1*numElementsLocal + i2*numElementsLocal*NE1;  
                      facePerm = face0;  
           
                      out->FaceElements->Id[k]=idCount++;  
                      out->FaceElements->Tag[k]=1;  
                      out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
                      out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;  
   
                      for( j=0; j<NUMNODES; j++ )  
                         out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
                  }  
              }  
              totalNECount+=NE1*NE2;  
              faceNECount+=NE1*NE2;  
      }  
      /* **  elements on boundary 002 (x1=1): */  
          if( domRight ) {  
              #pragma omp parallel for private(i1,i2,k)  
              for (i2=0;i2<NE2;i2++) {  
                  for (i1=0;i1<NE1;i1++) {  
                      k=i1+NE1*i2+faceNECount;  
                      kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;  
                      facePerm = face1;  
           
                      out->FaceElements->Id[k]=idCount++;  
                      out->FaceElements->Tag[k]=2;  
              out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
                      out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;  
   
                      for( j=0; j<NUMNODES; j++ )  
                         out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
                  }  
              }  
        totalNECount+=NE1*NE2;  
        faceNECount+=NE1*NE2;  
      }  
   }  
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x2=0): */  
     
      #pragma omp parallel for private(i0,i2,k)  
      for (i2=0;i2<NE2;i2++) {  
        for (i0=0;i0<numElementsLocal;i0++) {  
          k=i0+numElementsLocal*i2+faceNECount;  
                  kk=i0+numElementsLocal*NE1*i2;  
                  facePerm = face4;  
     
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Tag[k]=10;  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
          out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;  
   
                  for( j=0; j<NUMNODES; j++ )  
                     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
        }  
      }  
        if( boundaryLeft ){  
             for( i2=0; i2<NE2; i2++ )  
                 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;  
             if( periodicLocal[0] )  
                 for( i2=0; i2<NE2; i2++ )  
                     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;  
          }  
        if( boundaryRight )  
             for( i2=0; i2<NE2; i2++ )  
                 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
      totalNECount+=numElementsLocal*NE2;  
      faceNECount+=numElementsLocal*NE2;  
     
      /* **  elements on boundary 020 (x2=1): */  
     
      #pragma omp parallel for private(i0,i2,k)  
      for (i2=0;i2<NE2;i2++) {  
        for (i0=0;i0<numElementsLocal;i0++) {  
          k=i0+numElementsLocal*i2+faceNECount;  
                  kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;  
                  facePerm = face5;  
     
          out->FaceElements->Tag[k]=20;  
          out->FaceElements->Id[k]=idCount++;  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
          out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;  
   
                  for( j=0; j<NUMNODES; j++ )  
                     out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
        }  
      }  
        if( boundaryLeft ){  
             for( i2=0; i2<NE2; i2++ )  
                 out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;  
             if( periodicLocal[0] )  
                 for( i2=0; i2<NE2; i2++ )  
                     out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;  
          }  
        if( boundaryRight )  
             for( i2=0; i2<NE2; i2++ )  
                 out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
      totalNECount+=numElementsLocal*NE2;  
      faceNECount+=numElementsLocal*NE2;  
   }  
     out->FaceElements->elementDistribution->numInternal = faceNECount;  
       
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=23;  
     out->FaceElements->numElements=faceNECount;  
       
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
   
   /* setup distribution info for other elements */  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     /* reorder the degrees of freedom */  
     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );  
   
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
   if( !Finley_MPI_noError(mpi_info) )  
   {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out);  
   if( !Finley_MPI_noError(mpi_info) )  
   {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
   
447    /* free up memory */    /* free up memory */
448    Paso_MPIInfo_dealloc( mpi_info );    Paso_MPIInfo_free( mpi_info );  
   
449    #ifdef Finley_TRACE    #ifdef Finley_TRACE
450    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
451    #endif    #endif
   if (Finley_noError()) {  
        if ( ! Finley_Mesh_isPrepared(out) ) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   }  
   return out;  
 }  
 #endif  
452    
453      return out;
454    }

Legend:
Removed from v.1222  
changed lines
  Added in v.1223

  ViewVC Help
Powered by ViewVC 1.1.26