/[escript]/branches/domexper/dudley/src/Mesh_hex20.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Mesh_hex20.c

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

revision 787 by bcumming, Wed Jul 26 01:46:45 2006 UTC revision 2748 by gross, Tue Nov 17 07:32:59 2009 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2009 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 21  Line 23 
23    
24  /**************************************************************/  /**************************************************************/
25    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$  
   
 /**************************************************************/  
   
26  #include "RectangularMesh.h"  #include "RectangularMesh.h"
27    
28  #ifdef PASO_MPI  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,
29  /* get the number of nodes/elements for domain with rank=rank, of size processors                                            double* Length,
30     where n is the total number of nodes/elements in the global domain */                                            bool_t* periodic,
31  static index_t domain_MODdim( index_t rank, index_t size, index_t n )                                            index_t order,
32  {                                            index_t reduced_order,
33    rank = size-rank-1;                                            bool_t useElementsOnFace,
34                                              bool_t useFullElementOrder,
35    if( rank < n%size )                                            bool_t useMacroElements,
36      return (index_t)floor(n/size)+1;                                            bool_t optimize)
   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, dim_t *DOFBoundary )  
37  {  {
38    index_t i0;    #define N_PER_E 2
39    dim_t numNodesGlobal = numElementsGlobal+1;    #define DIM 3
40      dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0, Nstride2=0, local_NE0, local_NE1, local_NE2;
41    (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;    dim_t totalNECount,faceNECount,NDOF0=0, NDOF1=0, NDOF2=0, NFaceElements=0, local_N0=0, local_N1=0, local_N2=0, NN;
42    if( rank<size-1 ) // add on node for right hand boundary    index_t node0, myRank, e_offset0, e_offset1, e_offset2, offset0=0, offset1=0, offset2=0, global_i0, global_i1, global_i2;
     (*numNodesLocal) += 1;  
   
   numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;  
   periodicLocal[0] = periodicLocal[1] = FALSE;  
   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]--;  
     }  
   }  
   nodesExternal[0]*=2;  
   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 );  
   firstNode[0] *= 2;  
     
   numDOFLocal[0] = numNodesLocal[0];  
   if( periodicLocal[0] )  
     numDOFLocal[0]--;  
     
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
   
   if( size>1 )  
   {  
     DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;  
     DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;  
   }  
   else  
   {  
     DOFBoundary[0] = DOFBoundary[1] = 0;  
   }  
 }  
 #endif  
 /**************************************************************/  
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Hex20_singleCPU(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace, Paso_MPIInfo *mpi_info) {  
 #else  
 Finley_Mesh* Finley_RectangularMesh_Hex20(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;  
43    Finley_Mesh* out;    Finley_Mesh* out;
44      Paso_MPIInfo *mpi_info = NULL;
45      Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
46    char name[50];    char name[50];
47      bool_t generateAllNodes= useFullElementOrder || useMacroElements;
48      #ifdef Finley_TRACE
49    double time0=Finley_timer();    double time0=Finley_timer();
50      #endif
51    
52      /* get MPI information */
53      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
54      if (! Finley_noError()) {
55            return NULL;
56      }
57      myRank=mpi_info->rank;
58    
59      /* set up the global dimensions of the mesh */
60    
61    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
62    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
63    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
64    N0=2*NE0+1;    N0=N_PER_E*NE0+1;
65    N1=2*NE1+1;    N1=N_PER_E*NE1+1;
66    N2=2*NE2+1;    N2=N_PER_E*NE2+1;
67    
68    if (N0<=MIN(N1,N2)) {    /*  allocate mesh: */  
69       if (N1 <= N2) {    sprintf(name,"Brick %d x %d x %d mesh",N0,N1,N2);
70          M0=1;    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
71          M1=N0;    if (! Finley_noError()) {
72          M2=N0*N1;        Paso_MPIInfo_free( mpi_info );
73       } else {        return NULL;
74          M0=1;    }
75          M2=N0;  
76          M1=N0*N2;    if (generateAllNodes) {
77       }       /* Finley_setError(SYSTEM_ERROR,"full element order for Hex elements is not supported yet."); */
78    } else if (N1<=MIN(N2,N0)) {       if (useMacroElements) {
79       if (N2 <= N0) {            refElements= Finley_ReferenceElementSet_alloc(Hex27Macro,out->order,out->reduced_order);
         M1=1;  
         M2=N1;  
         M0=N2*N1;  
80       } else {       } else {
81          M1=1;            refElements=Finley_ReferenceElementSet_alloc(Hex27, out->order,out->reduced_order);
         M0=N1;  
         M2=N1*N0;  
82       }       }
83    } else {       if (useElementsOnFace) {
84       if (N0 <= N1) {           Finley_setError(SYSTEM_ERROR,"rich elements for Hex27 elements is not supported yet.");
         M2=1;  
         M0=N2;  
         M1=N2*N0;  
85       } else {       } else {
86          M2=1;           if (useMacroElements) {
87          M1=N2;               refFaceElements=Finley_ReferenceElementSet_alloc(Rec9Macro, out->order, out->reduced_order);
88          M0=N1*N2;           } else {
89                 refFaceElements=Finley_ReferenceElementSet_alloc(Rec9, out->order, out->reduced_order);
90             }
91             refContactElements=Finley_ReferenceElementSet_alloc(Rec9_Contact, out->order, out->reduced_order);
92       }       }
   }  
93    
94    NFaceElements=0;    } else  {
95    if (!periodic[0]) {       refElements= Finley_ReferenceElementSet_alloc(Hex20,out->order,out->reduced_order);
96        NDOF0=N0;       if (useElementsOnFace) {
97        NFaceElements+=2*NE1*NE2;           refFaceElements = Finley_ReferenceElementSet_alloc(Hex20Face ,out->order,out->reduced_order);
98    } else {           refContactElements=Finley_ReferenceElementSet_alloc(Hex20Face_Contact, out->order, out->reduced_order);
       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;  
   } else {  
       NDOF2=N2-1;  
   }  
     
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);  
 #ifdef PASO_MPI  
   out=Finley_Mesh_alloc(name,3,order,mpi_info);  
   
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Hex20,out->order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
99    
100           } else {
101    /*  allocate tables: */           refFaceElements = Finley_ReferenceElementSet_alloc(Rec8 ,out->order,out->reduced_order);
102    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);           refContactElements=Finley_ReferenceElementSet_alloc(Rec8_Contact, out->order, out->reduced_order);
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
 #else  
   out=Finley_Mesh_alloc(name,3,order);  
   
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Hex20,out->order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   
     /* create 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  
       }  
     }  
   }  
103    
     
   /* 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;  
        }  
104       }       }
105    }    }
106    if (!periodic[1]) {    refPoints=Finley_ReferenceElementSet_alloc(Point1, out->order, out->reduced_order);
107      for (i2=0;i2<N2;i2++) {  
108        for (i0=0;i0<N0;i0++) {    if ( Finley_noError()) {
          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=2*i0+2*i1*N0+2*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  
   
         out->Elements->Nodes[INDEX2(0,k,20)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;  
         out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;  
         out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;  
         out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;  
         out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;  
         out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;  
         out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;  
         out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;  
         out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;  
         out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;  
         out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;  
         out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;  
         out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;  
         out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;  
         out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;  
         out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;  
         out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;  
         out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;  
         out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;  
       }  
     }  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);  
     
   /*   face elements: */  
109        
110    if  (useElementsOnFace) {        Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
111       NUMNODES=20;        Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
112    } else {        Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
113       NUMNODES=8;        Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
114    
115          /* work out the largest dimension */
116          if (N2==MAX3(N0,N1,N2)) {
117              Nstride0=1;
118              Nstride1=N0;
119              Nstride2=N0*N1;
120              local_NE0=NE0;
121              e_offset0=0;
122              local_NE1=NE1;
123              e_offset1=0;
124              Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
125          } else if (N1==MAX3(N0,N1,N2)) {
126              Nstride0=N2;
127              Nstride1=N0*N2;
128              Nstride2=1;
129              local_NE0=NE0;
130              e_offset0=0;
131              Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
132              local_NE2=NE2;
133              e_offset2=0;
134          } else {
135              Nstride0=N1*N2;
136              Nstride1=1;
137              Nstride2=N1;
138              Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
139              local_NE1=NE1;
140              e_offset1=0;
141              local_NE2=NE2;
142              e_offset2=0;
143          }
144          offset0=e_offset0*N_PER_E;
145          offset1=e_offset1*N_PER_E;
146          offset2=e_offset2*N_PER_E;
147          local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
148          local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
149          local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
150    
151          /* get the number of surface elements */
152    
153          NFaceElements=0;  
154          if (!periodic[2] && (local_NE2>0) ) {
155              NDOF2=N2;
156              if (offset2==0) NFaceElements+=local_NE1*local_NE0;
157              if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
158          } else {
159              NDOF2=N2-1;
160          }
161    
162          if (!periodic[0] && (local_NE0>0) ) {
163              NDOF0=N0;
164              if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
165              if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
166          } else {
167              NDOF0=N0-1;
168          }
169          if (!periodic[1] && (local_NE1>0) ) {
170              NDOF1=N1;
171              if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
172              if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
173          } else {
174              NDOF1=N1-1;
175          }
176    }    }
177    totalNECount=NE0*NE1*NE2;    /*  allocate tables: */
178    faceNECount=0;    if (Finley_noError()) {
179            Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
180    /*   these are the quadrilateral elements on boundary 1 (x3=0): */        Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
181            Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
   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=2*i0+2*i1*N0;  
     
         out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;  
         out->FaceElements->Tag[k]=100;  
         out->FaceElements->Color[k]=(i0%2)+2*(i1%2);  
 #ifdef PASO_MPI  
         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
     
         if  (useElementsOnFace) {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;  
            out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;  
            out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;  
            out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;  
            out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;  
         } else {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
         }  
       }  
     }  
     totalNECount+=NE1*NE0;  
     faceNECount+=NE1*NE0;  
     
     /* **  elements on boundary 200 (x3=1) */  
     #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=2*i0+2*i1*N0+2*N0*N1*(NE2-1);  
     
         out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;  
         out->FaceElements->Tag[k]=200;  
         out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;  
 #ifdef PASO_MPI  
         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
         if  (useElementsOnFace) {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;  
            out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;  
   
            out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;  
            out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;  
            out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;  
   
         } else {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;  
         }  
       }  
     }  
     totalNECount+=NE1*NE0;  
     faceNECount+=NE1*NE0;  
182    }    }
183    if (!periodic[0]) {    if (Finley_noError()) {
      /* **  elements on boundary 001 (x1=0): */  
     
      #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=2*i1*N0+2*N0*N1*i2;  
     
          out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;  
          out->FaceElements->Tag[k]=1;  
          out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;  
 #ifdef PASO_MPI  
         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
   
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;  
   
             out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;  
   
             out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;  
             out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
             out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;  
   
             out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
             out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;  
184    
185           } else {        /* create nodes */
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
          }  
        }  
      }  
      totalNECount+=NE1*NE2;  
      faceNECount+=NE1*NE2;  
     
      /* **  elements on boundary 002 (x1=1): */  
     
      #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=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;  
186        
187           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;            #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
188           out->FaceElements->Tag[k]=2;       for (i2=0;i2<local_N2;i2++) {
189           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;         for (i1=0;i1<local_N1;i1++) {
190  #ifdef PASO_MPI           for (i0=0;i0<local_N0;i0++) {
191          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             k=i0+local_N0*i1+local_N0*local_N1*i2;
192  #endif             global_i0=i0+offset0;
193               global_i1=i1+offset1;
194           if  (useElementsOnFace) {             global_i2=i2+offset2;
195              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
196              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
197              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;             out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
198              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
199               out->Nodes->Tag[k]=0;
200              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
201              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;                                                 +Nstride1*(global_i1%NDOF1)
202              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;                                                 +Nstride2*(global_i2%NDOF2);
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;  
   
             out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;  
             out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
             out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;  
   
             out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;  
             out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
             out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;  
   
             out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;  
   
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;  
203           }           }
     
204         }         }
205       }       }
206       totalNECount+=NE1*NE2;       /*   set the elements: */
207       faceNECount+=NE1*NE2;       NN=out->Elements->numNodes;
208    }       #pragma omp parallel for private(i0,i1,i2,k,node0)
209    if (!periodic[1]) {       for (i2=0;i2<local_NE2;i2++) {
210       /* **  elements on boundary 010 (x2=0): */         for (i1=0;i1<local_NE1;i1++) {
211               for (i0=0;i0<local_NE0;i0++) {
212       #pragma omp parallel for private(i0,i2,k,node0)            
213       for (i2=0;i2<NE2;i2++) {             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
214         for (i0=0;i0<NE0;i0++) {             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
          k=i0+NE0*i2+faceNECount;  
          node0=2*i0+2*N0*N1*i2;  
215        
216           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
217           out->FaceElements->Tag[k]=10;             out->Elements->Tag[k]=0;
218           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;             out->Elements->Owner[k]=myRank;
219  #ifdef PASO_MPI  
220          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             out->Elements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
221  #endif             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                      2*Nstride0;
222               out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           2*Nstride1+2*Nstride0;
223           if  (useElementsOnFace) {             out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride1;
224              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2                      ;
225              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;             out->Elements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2           +2*Nstride0;
226              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
227              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
228               out->Elements->Nodes[INDEX2(8,k,NN)] =node0+                      1*Nstride0;
229              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;             out->Elements->Nodes[INDEX2(9,k,NN)] =node0+           1*Nstride1+2*Nstride0;
230              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;             out->Elements->Nodes[INDEX2(10,k,NN)]=node0+           2*Nstride1+1*Nstride0;
231              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;             out->Elements->Nodes[INDEX2(11,k,NN)]=node0+           1*Nstride1           ;
232              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;             out->Elements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2                      ;
233               out->Elements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2           +2*Nstride0;
234              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;             out->Elements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
235              out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;             out->Elements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
236              out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;             out->Elements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2           +1*Nstride0;
237              out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;             out->Elements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
238               out->Elements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
239              out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;             out->Elements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
240              out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;             if (generateAllNodes) {
241              out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;                out->Elements->Nodes[INDEX2(20,k,NN)]=node0+           1*Nstride1+1*Nstride0;
242              out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;                out->Elements->Nodes[INDEX2(21,k,NN)]=node0+1*Nstride2           +1*Nstride0;
243                  out->Elements->Nodes[INDEX2(22,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
244              out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;                out->Elements->Nodes[INDEX2(23,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
245              out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;                out->Elements->Nodes[INDEX2(24,k,NN)]=node0+1*Nstride2+1*Nstride1           ;
246              out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;                out->Elements->Nodes[INDEX2(25,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
247              out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;                out->Elements->Nodes[INDEX2(26,k,NN)]=node0+1*Nstride2+1*Nstride1+1*Nstride0;        
248               }
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;  
249           }           }
250         }         }
251       }       }
252       totalNECount+=NE0*NE2;       /* face elements */
253       faceNECount+=NE0*NE2;       NN=out->FaceElements->numNodes;
254           totalNECount=NE0*NE1*NE2;
255       /* **  elements on boundary 020 (x2=1): */       faceNECount=0;
256         /*   these are the quadrilateral elements on boundary 1 (x3=0): */
257         if (!periodic[2] && (local_NE2>0)) {
258           /* **  elements on boundary 100 (x3=0): */
259           if (offset2==0) {
260              #pragma omp parallel for private(i0,i1,k,node0)
261              for (i1=0;i1<local_NE1;i1++) {
262                for (i0=0;i0<local_NE0;i0++) {
263              
264                  k=i0+local_NE0*i1+faceNECount;
265                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
266            
267       #pragma omp parallel for private(i0,i2,k,node0)                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
268       for (i2=0;i2<NE2;i2++) {                out->FaceElements->Tag[k]=100;
269         for (i0=0;i0<NE0;i0++) {                out->FaceElements->Owner[k]=myRank;
270           k=i0+NE0*i2+faceNECount;          
271           node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;                if  (useElementsOnFace) {
272                       out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
273           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0           +2*Nstride1           ;
274           out->FaceElements->Tag[k]=20;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0           +2*Nstride1+2*Nstride0;
275           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride0           ;
276  #ifdef PASO_MPI                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2                      ;
277          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
278  #endif                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
279                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2           +2*Nstride0;
280           if  (useElementsOnFace) {                   out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+           1*Nstride1           ;
281              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;                   out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+           2*Nstride1+1*Nstride0;
282              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;                   out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+           1*Nstride1+2*Nstride0;
283              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;                   out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+                      1*Nstride0;
284              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;                   out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2                      ;
285                     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
286              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
287              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;                   out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2           +2*Nstride0;
288              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;                   out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2+1*Nstride1;
289              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;                   out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
290                     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
291              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;                   out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2           +1*Nstride0;
292              out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;                } else {
293              out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
294              out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           2*Nstride1           ;
295                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           2*Nstride1+2*Nstride0;
296              out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                      2*Nstride0;
297              out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+           1*Nstride1           ;
298              out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+           2*Nstride1+1*Nstride0;
299              out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           1*Nstride1+2*Nstride0;
300                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+                      1*Nstride0;
301              out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;                   if (generateAllNodes){
302              out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+           1*Nstride1+1*Nstride0;
303              out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;                   }
304              out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;                }
305           } else {              }
306              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;            }
307              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;            faceNECount+=local_NE1*local_NE0;
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;  
          }  
308         }         }
309       }         totalNECount+=NE1*NE0;
310       totalNECount+=NE0*NE2;         /* **  elements on boundary 200 (x3=1) */
311       faceNECount+=NE0*NE2;         if (local_NE2+e_offset2 == NE2) {
312    }            #pragma omp parallel for private(i0,i1,k,node0)
313    out->FaceElements->minColor=0;            for (i1=0;i1<local_NE1;i1++) {
314    out->FaceElements->maxColor=24;              for (i0=0;i0<local_NE0;i0++) {
315          
316  #ifdef PASO_MPI                k=i0+local_NE0*i1+faceNECount;
317      Finley_ElementFile_setDomainFlags( out->Elements );                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
318      Finley_ElementFile_setDomainFlags( out->FaceElements );          
319      Finley_ElementFile_setDomainFlags( out->ContactElements );                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
320      Finley_ElementFile_setDomainFlags( out->Points );                out->FaceElements->Tag[k]=200;
321                  out->FaceElements->Owner[k]=myRank;
322      /* reorder the degrees of freedom */                if  (useElementsOnFace) {
323      Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2                      ;
324  #endif                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2+           2*Nstride0;
325                         out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
326    /*   condense the nodes: */                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
327    Finley_Mesh_resolveNodeIds(out);        
328                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
329    /* prepare mesh for further calculatuions:*/                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride0                      ;
330    Finley_Mesh_prepare(out) ;                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           2*Nstride1+2*Nstride0;
331                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           2*Nstride1;
332    #ifdef Finley_TRACE        
333    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);                   out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+           1*Nstride0;
334    #endif                   out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
335                     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
336    if (! Finley_noError()) {                   out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
337        Finley_Mesh_dealloc(out);        
338        return NULL;                   out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
339    }                   out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2           +2*Nstride0;
340    return out;                   out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
341  }                   out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
342  #ifdef PASO_MPI        
343  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {                   out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+                      1*Nstride0;
344    dim_t N0,N0t,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF0t,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;                   out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+           1*Nstride1+2*Nstride0;
345    dim_t kk,iI, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];                   out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+           2*Nstride1+1*Nstride0;
346      bool_t dom_left, dom_right, dom_internal;                   out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+           1*Nstride1           ;
347    index_t firstNode=0, DOFcount=0, node0, node1, node2, idCount;        
348    index_t targetDomain=-1, firstNodeConstruct, j;                } else {
349    bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2                      ;
350      index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2           +2*Nstride0;
351    Finley_Mesh* out;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
352    char name[50];                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
353    double time0=Finley_timer();                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2           +1*Nstride0;
354    Paso_MPIInfo *mpi_info = NULL;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
355                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+1*Nstride0;
356    NE0=MAX(1,numElements[0]);                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
357    NE1=MAX(1,numElements[1]);                   if (generateAllNodes){
358    NE2=MAX(1,numElements[2]);                   out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+1*Nstride1+1*Nstride0;
359    N0=2*NE0+1;                   }
360    N1=2*NE1+1;                }
361    N2=2*NE2+1;              }
362              }
363      index_t face0[] = {0,4,7,3,1,5,6,2,12,19,15,11,8,16,18,10,13,17,14,9};            faceNECount+=local_NE1*local_NE0;
     index_t face1[] = {1,2,6,5,0,3,7,4,9,14,17,13,8,10,18,16,11,15,19,12};  
     index_t face2[] = {0,1,5,4,3,2,6,7,8,13,16,12,11,9,17,19,14,18,15,10};  
     index_t face3[] = {3,7,6,2,0,4,5,1,15,18,14,10,11,19,17,9,12,16,13,8};  
     index_t face4[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,19,18,17,16};  
     index_t face5[] = {4,5,6,7,0,1,2,3,16,17,18,19,12,13,14,15,8,9,10,11};  
       
     index_t face0R[] = {0,4,7,3,12,19,15,11};  
     index_t face1R[] = {1,2,6,5,9,14,17,13};  
     index_t face2R[] = {0,1,5,4,8,13,16,12};  
     index_t face3R[] = {3,7,6,2,15,18,14,10};  
     index_t face4R[] = {0,3,2,1,11,10,9,8};  
     index_t face5R[] = {4,5,6,7,16,17,18,19};  
       
   /* 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_Hex20_singleCPU( numElements, Length, periodic, 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, DOFBoundary );    
   
   NFaceElements=0;  
   if (!periodic[0]) {  
       NDOF0=N0;  
       NFaceElements+=(domRight+domLeft)*NE1*NE2;  
   } else {  
       NDOF0=N0-1;  
   }  
   if (!periodic[1]) {  
       NDOF1=N1;  
       NFaceElements+=2*numElementsLocal*NE2;  
   } else {  
       NDOF1=N1-1;  
   }  
   if (!periodic[2]) {  
       NDOF2=N2;  
       NFaceElements+=2*numElementsLocal*NE1;  
   } else {  
       NDOF2=N2-1;  
   }  
     
     boundaryLeft = !domLeft || periodicLocal[0];  
     boundaryRight = !domRight || periodicLocal[1];  
     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;  
     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;  
     firstNodeConstruct = firstNode - 2*boundaryLeft;  
     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;  
   
   /*  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(Hex20,out->order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8_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,N0t*N1*N2);  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*3, 0 );  
   Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1*NE2);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
   
     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_EXTERNAL;  
                 out->Nodes->Dom[N1*N0t*i2+N0t*i1+2] = 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;  
                 out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-3] = 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[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;  
364         }         }
365           totalNECount+=NE1*NE0;
366       }       }
367    }       if (!periodic[0] && (local_NE0>0)) {
368    if (!periodic[1]) {          /* **  elements on boundary 001 (x1=0): */
369      for (i2=0;i2<N2;i2++) {      
370        for (i0=0;i0<N0t;i0++) {          if (e_offset0 == 0) {
371           out->Nodes->Tag[i0 + i2*N1*N0t]+=10;             #pragma omp parallel for private(i1,i2,k,node0)
372           out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;             for (i2=0;i2<local_NE2;i2++) {
373        }               for (i1=0;i1<local_NE1;i1++) {
374      }        
375    }                 k=i1+local_NE1*i2+faceNECount;
376    if (!periodic[0] && !domInternal ) {                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
377      for (i2=0;i2<N2;i2++) {                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
378        for (i1=0;i1<N1;i1++) {                 out->FaceElements->Tag[k]=1;
379          if( domLeft )                 out->FaceElements->Owner[k]=myRank;
380            out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;        
381          if( domRight )                 if  (useElementsOnFace) {
382            out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
383        }                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2                      ;
384      }                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
385    }                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride1                      ;
386          
387      /* form the boudary communication information */                    out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride0                      ;
388      forwardDOF  = MEMALLOC(NDOF1*NDOF2*2,index_t);                    out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride0           ;
389      backwardDOF = MEMALLOC(NDOF1*NDOF2*2,index_t);                    out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
390      if( !(mpi_info->size==2 && periodicLocal[0])){                    out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride1+2*Nstride0           ;
391          if( boundaryLeft  ) {        
392              targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;                    out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2                      ;
393              for( iI=0, i2=0; i2<NDOF2; i2++ ){                    out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
394                  for( i1=0; i1<NDOF1; i1++, iI+=2 ){                    out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
395                      forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];                    out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+           1*Nstride1           ;
396                      backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];        
397                      backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];                    out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+                      1*Nstride0;
398                  }                    out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2           +1*Nstride0;
399              }                    out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
400              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );                    out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride1+           1*Nstride0;
401              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );        
402          }                    out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2+           2*Nstride0;
403          if( boundaryRight ) {                    out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
404              targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;                    out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
405              for( iI=0,i2=0; i2<NDOF2; i2++ ){                    out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride1+           2*Nstride0;
406                  for( i1=0; i1<NDOF1; i1++, iI+=2 ){        
407                      forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];                 } else {
408                      forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];                    out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
409                      backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];                    out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2                      ;
410                  }                    out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
411              }                    out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride1           ;
412              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );                    out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+1*Nstride2                      ;
413              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );                    out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
414          }                    out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+1*Nstride2+2*Nstride1           ;
415      } else{                    out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           1*Nstride1           ;
416          /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */                   if (generateAllNodes){
417          targetDomain = 1;                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1           ;
418                             }
419          for( iI=0,i2=0; i2<NDOF2; i2++ ){                 }
420              for( i1=0; i1<NDOF1; i1++, iI+=2 ){               }
421                  forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];             }
422                  forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];             faceNECount+=local_NE1*local_NE2;
423                  backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];          }
424              }          totalNECount+=NE1*NE2;
425          }      
426          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );          /* **  elements on boundary 002 (x1=1): */
427          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );          if (local_NE0+e_offset0 == NE0) {
428               #pragma omp parallel for private(i1,i2,k,node0)
429          for( iI=0, i2=0; i2<NDOF2; i2++ ){             for (i2=0;i2<local_NE2;i2++) {
430              for( i1=0; i1<NDOF1; i1++, iI+=2 ){               for (i1=0;i1<local_NE1;i1++) {
431                  forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];                 k=i1+local_NE1*i2+faceNECount;
432                  backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];                 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
433                  backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
434              }                 out->FaceElements->Tag[k]=2;
435          }                 out->FaceElements->Owner[k]=myRank;
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );  
     }  
     MEMFREE( forwardDOF );  
     MEMFREE( backwardDOF );  
   /*   set the 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) ? 2*(i1*N0t + i2*N1*N0t) : 2*(i1*N0t + i2*N1*N0t + i0) + periodicLocal[0];  
   
                 out->Elements->Id[k]=((firstNodeConstruct/2+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,20)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;  
         out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0t+2;  
         out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0t;  
         out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0t*N1;  
         out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0t*N1+2;  
         out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0t*N1+2*N0t+2;  
         out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0t*N1+2*N0t;  
         out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;  
         out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0t+2;  
         out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0t+1;  
         out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0t;  
         out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0t*N1;  
         out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0t*N1+2;  
         out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0t*N1+2*N0t+2;  
         out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0t*N1+2*N0t;  
         out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0t*N1+1;  
         out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0t*N1+N0t+2;  
         out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0t*N1+2*N0t+1;  
         out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*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;  
   
     out->Elements->numElements = numElementsLocal*NE1*NE2;  
     Finley_ElementFile_setDomainFlags( out->Elements );  
     
   /*   face elements: */  
     
   if  (useElementsOnFace) {  
      NUMNODES=20;  
   } else {  
      NUMNODES=8;  
   }  
   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[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 =!useElementsOnFace ? face0R :  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,20)];  
                  }  
              }  
              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 =!useElementsOnFace ? face1R :  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,20)];  
                  }  
              }  
        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 =!useElementsOnFace ? face2R :  face2;  
436        
437           out->FaceElements->Id[k]=idCount++;                 if  (useElementsOnFace) {
438           out->FaceElements->Tag[k]=10;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      2*Nstride0;
439           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           2*Nstride1+2*Nstride0;
440           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
441                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+           2*Nstride0;
442                   for( j=0; j<NUMNODES; j++ )        
443                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
444         }                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           2*Nstride1           ;
445                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
446                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2                      ;
447          
448                      out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+           1*Nstride1+2*Nstride0;
449                      out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
450                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
451                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2+           2*Nstride0;
452          
453                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+                      1*Nstride0;
454                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+           2*Nstride1+1*Nstride0;
455                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
456                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+           1*Nstride0;
457          
458                      out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+           1*Nstride1           ;
459                      out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
460                      out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
461                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2                      ;
462          
463                   } else {
464                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +2*Nstride0;
465                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           2*Nstride1+2*Nstride0;
466                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
467                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+           2*Nstride0;
468                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           1*Nstride1+2*Nstride0;
469                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
470                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
471                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2           +2*Nstride0;
472                     if (generateAllNodes){
473                        out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1+2*Nstride0;
474                     }
475                   }
476            
477                 }
478               }
479               faceNECount+=local_NE1*local_NE2;
480             }
481             totalNECount+=NE1*NE2;
482       }       }
483         if( boundaryLeft ){       if (!periodic[1] && (local_NE1>0)) {
484              for( i2=0; i2<NE2; i2++ )          /* **  elements on boundary 010 (x2=0): */
485                  out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;          if (e_offset1 == 0) {
486              if( periodicLocal[0] )             #pragma omp parallel for private(i0,i2,k,node0)
487                  for( i2=0; i2<NE2; i2++ )             for (i2=0;i2<local_NE2;i2++) {
488                      out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;               for (i0=0;i0<local_NE0;i0++) {
489           }                 k=i0+local_NE0*i2+faceNECount;
490         if( boundaryRight )                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
491              for( i2=0; i2<NE2; i2++ )          
492                  out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
493       totalNECount+=numElementsLocal*NE2;                 out->FaceElements->Tag[k]=10;
494       faceNECount+=numElementsLocal*NE2;                 out->FaceElements->Owner[k]=myRank;
495                     if  (useElementsOnFace) {
496       /* **  elements on boundary 020 (x2=1): */                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
497                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      2*Nstride0;
498       #pragma omp parallel for private(i0,i2,k)                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2           +2*Nstride0;
499       for (i2=0;i2<NE2;i2++) {                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2                      ;
500         for (i0=0;i0<numElementsLocal;i0++) {        
501           k=i0+numElementsLocal*i2+faceNECount;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           2*Nstride1           ;
502                   kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+           2*Nstride0;
503                   facePerm =!useElementsOnFace ? face3R :  face3;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
504                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
505          
506                      out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+                      1*Nstride0;
507                      out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+           2*Nstride0;
508                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+           1*Nstride0;
509                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2                      ;
510          
511                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+           1*Nstride1           ;
512                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+           1*Nstride1+2*Nstride0;
513                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
514                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
515        
516           out->FaceElements->Tag[k]=20;                    out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+           2*Nstride1+1*Nstride0;
517           out->FaceElements->Id[k]=idCount++;                    out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
518           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
519           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;                    out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
520          
521                   for( j=0; j<NUMNODES; j++ )                 } else {
522                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
523         }                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      2*Nstride0;
524       }                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+           2*Nstride0;
525         if( boundaryLeft ){                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2                      ;
526              for( i2=0; i2<NE2; i2++ )                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+                      1*Nstride0;
527                  out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+           2*Nstride0;
528              if( periodicLocal[0] )                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+           1*Nstride0;
529                  for( i2=0; i2<NE2; i2++ )                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2                      ;
530                      out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;                   if (generateAllNodes){
531           }                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+         1*Nstride0;
532         if( boundaryRight )                   }
533              for( i2=0; i2<NE2; i2++ )                 }
534                  out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;               }
535       totalNECount+=numElementsLocal*NE2;             }
536       faceNECount+=numElementsLocal*NE2;             faceNECount+=local_NE0*local_NE2;
537    }          }
538    if (!periodic[2]) {          totalNECount+=NE0*NE2;
539       /*  elements on boundary 100 (x3=0): */          /* **  elements on boundary 020 (x2=1): */
540              if (local_NE1+e_offset1 == NE1) {
541       #pragma omp parallel for private(i0,i1,k)             #pragma omp parallel for private(i0,i2,k,node0)
542       for (i1=0;i1<NE1;i1++) {             for (i2=0;i2<local_NE2;i2++) {
543         for (i0=0; i0<numElementsLocal; i0++) {               for (i0=0;i0<local_NE0;i0++) {
544           k=i0+numElementsLocal*i1+faceNECount;                 k=i0+local_NE0*i2+faceNECount;
545                   kk=i0 + i1*numElementsLocal;                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
                  facePerm =!useElementsOnFace ? face4R :  face4;  
546        
547           out->FaceElements->Id[k]=idCount++;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
548           out->FaceElements->Tag[k]=100;                 out->FaceElements->Tag[k]=20;
549           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                 out->FaceElements->Owner[k]=myRank;
550           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);        
551                   if  (useElementsOnFace) {
552                   for( j=0; j<NUMNODES; j++ )                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           2*Nstride1           ;
553                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
554         }                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
555                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0           ;
556          
557                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
558                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2                      ;
559                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+           2*Nstride0;
560                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      2*Nstride0;
561          
562                      out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
563                      out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
564                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
565                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+           2*Nstride1+1*Nstride0;
566          
567                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+           1*Nstride1           ;
568                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
569                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
570                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+           1*Nstride1+2*Nstride0;
571          
572                      out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2                      ;
573                      out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2           +1*Nstride0;
574                      out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2           +2*Nstride0;
575                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+                      1*Nstride0;
576                   } else {
577                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           2*Nstride1           ;
578                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
579                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
580                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           2*Nstride1+2*Nstride0;
581                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
582                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
583                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
584                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+           2*Nstride1+1*Nstride0;
585                     if (generateAllNodes){
586                        out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+2*Nstride1+1*Nstride0;
587                     }
588                   }
589                 }
590               }
591               faceNECount+=local_NE0*local_NE2;
592            }
593            totalNECount+=NE0*NE2;
594       }       }
595         if( boundaryLeft ){       if (Finley_noError()) {
596              for( i1=0; i1<NE1; i1++ )           /* add tag names */
597                  out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;           Finley_Mesh_addTagMap(out,"top", 200);
598              if( periodicLocal[0] )           Finley_Mesh_addTagMap(out,"bottom", 100);
599                  for( i1=0; i1<NE1; i1++ )           Finley_Mesh_addTagMap(out,"left", 1);
600                      out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;           Finley_Mesh_addTagMap(out,"right", 2);
601           }           Finley_Mesh_addTagMap(out,"front", 10);
602         if( boundaryRight )           Finley_Mesh_addTagMap(out,"back", 20);
             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 = !useElementsOnFace ? face5R : face5;  
   
          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,20)];  
        }  
603       }       }
        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;  
604    }    }
     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;  
   }  
   
605    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
606    Finley_Mesh_prepare(out);    if (Finley_noError()) {
607    if( !Finley_MPI_noError(mpi_info) )           Finley_Mesh_resolveNodeIds(out);
608    {    }
609      Paso_MPIInfo_dealloc( mpi_info );    if (Finley_noError()) {
610      Finley_Mesh_dealloc(out);           Finley_Mesh_prepare(out, optimize);
611      return NULL;    }
   }  
   
   /* free up memory */  
   Paso_MPIInfo_dealloc( mpi_info );  
   
     //print_mesh_statistics( out, FALSE );  
612    
613    #ifdef Finley_TRACE    if (!Finley_noError()) {
614    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);        Finley_Mesh_free(out);
615    #endif    }
616        /* free up memory */
617      Finley_ReferenceElementSet_dealloc(refPoints);
618      Finley_ReferenceElementSet_dealloc(refContactElements);
619      Finley_ReferenceElementSet_dealloc(refFaceElements);
620      Finley_ReferenceElementSet_dealloc(refElements);
621      Paso_MPIInfo_free( mpi_info );  
622    
623    return out;    return out;
624  }  }
 #endif  
   

Legend:
Removed from v.787  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26