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

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

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

revision 730 by bcumming, Mon May 15 04:03:49 2006 UTC revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC
# Line 29  Line 29 
29    
30  /**************************************************************/  /**************************************************************/
31    
32  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {  #ifdef PASO_MPI
33    /* get the number of nodes/elements for domain with rank=rank, of size processors
34       where n is the total number of nodes/elements in the global domain */
35    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36    {
37      rank = size-rank-1;
38    
39      if( rank < n%size )
40        return (index_t)floor(n/size)+1;
41      return (index_t)floor(n/size);
42    }
43    
44    
45    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46    /* A bit messy, but it only has to be done once... */
47    static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal )
48    {
49      index_t i0;
50      dim_t numNodesGlobal = numElementsGlobal+1;
51    
52      (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53      
54      numElementsLocal[0] = numNodesLocal[0]+1;
55      periodicLocal[0] = periodicLocal[1] = FALSE;
56      nodesExternal[0] = nodesExternal[1] = 1;
57      if( periodic )
58      {
59        if( size==1 )
60        {
61          numElementsLocal[0] = numElementsGlobal;
62          nodesExternal[0] = nodesExternal[1] = 0;
63          periodicLocal[0] = periodicLocal[1] = TRUE;
64        }
65        else
66        {
67          if( rank==0 )
68          {
69            periodicLocal[0] = TRUE;
70            numNodesLocal[0]++;
71          }
72          if( rank==(size-1) )
73          {      
74            periodicLocal[1] = TRUE;
75            numNodesLocal[0]--;
76            numElementsLocal[0]--;
77          }
78        }
79      }
80      else if( !periodic )
81      {
82        if( rank==0 ){
83          nodesExternal[0]--;
84          numElementsLocal[0]--;
85        }
86        if( rank==(size-1) )
87        {
88          nodesExternal[1]--;
89          numElementsLocal[0]--;
90        }
91      }
92      numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93      
94      numElementsInternal[0] = numElementsLocal[0];
95      if( (rank==0) && (rank==size-1) );
96        
97      else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98          numElementsInternal[0] -= 1;
99      else
100        numElementsInternal[0] -= 2;
101    
102      firstNode[0] = 0;
103      for( i0=0; i0<rank; i0++ )
104        firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105    
106      numDOFLocal[0] = numNodesLocal[0];
107      if( periodicLocal[0] )
108      {
109        numDOFLocal[0]--;
110      }
111      DOFExternal[0] = nodesExternal[0];
112      DOFExternal[1] = nodesExternal[1];
113    }
114    
115    void print_mesh_statistics( Finley_Mesh *out )
116    {
117      index_t i, j, N;
118      
119      printf( "\nNodes\n=====\n\n" );
120        printf( "\t%d internal DOF\n\t%d boundary DOF\n\t%d local DOF\n\t%d external DOF\n", out->Nodes->degreeOfFreedomDistribution->numInternal, out->Nodes->degreeOfFreedomDistribution->numBoundary, out->Nodes->degreeOfFreedomDistribution->numLocal, out->Nodes->degreeOfFreedomDistribution->numExternal);
121      for( i=0; i<out->Nodes->numNodes; i++ )
122        printf( "node %d\t: id %d   \tDOF %d   \t: tag %d  \t: coordinates [%3g, %3g, %3g]\n", i, out->Nodes->Id[i], out->Nodes->degreeOfFreedom[i], out->Nodes->Tag[i], out->Nodes->Coordinates[INDEX2(0,i,3)], out->Nodes->Coordinates[INDEX2(1,i,3)], out->Nodes->Coordinates[INDEX2(2,i,3)] );
123    
124      printf( "Elements\n========\n\n" );
125        printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->Elements->elementDistribution->numInternal, out->Elements->elementDistribution->numBoundary, out->Elements->elementDistribution->numLocal );
126        N = out->Elements->ReferenceElement->Type->numNodes;
127      for( i=0; i<out->Elements->numElements; i++ ){
128        printf( "element %d    \t: id %d  \t: nodes [ %3d, %3d, %3d, %3d, %3d, %3d, %3d, %3d ]", i, out->Elements->Id[i], out->Elements->Nodes[INDEX2(0,i,8)], out->Elements->Nodes[INDEX2(1,i,8)], out->Elements->Nodes[INDEX2(2,i,8)], out->Elements->Nodes[INDEX2(3,i,8)], out->Elements->Nodes[INDEX2(4,i,8)], out->Elements->Nodes[INDEX2(5,i,8)], out->Elements->Nodes[INDEX2(6,i,8)], out->Elements->Nodes[INDEX2(7,i,8)] );
129            printf( " DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(0,i,N)]]] );  
130            for( j=1; j<N; j++ )
131                printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(j,i,N)]]]  );  
132            printf( " ]\n" );  
133      }
134    
135        printf( "\nFace Elements\n==============\n\n" );
136        printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->FaceElements->elementDistribution->numInternal, out->FaceElements->elementDistribution->numBoundary, out->FaceElements->elementDistribution->numLocal );
137        N = out->FaceElements->ReferenceElement->Type->numNodes;
138      for( i=0; i<out->FaceElements->numElements; i++ ){
139            printf( "face element %d \t: id %d  \t: nodes [ %3d", i, out->FaceElements->Id[i], out->FaceElements->Nodes[INDEX2(0,i,N)] );
140            for( j=1; j<N; j++ )
141                printf( ", %3d", out->FaceElements->Nodes[INDEX2(j,i,N)]  );    
142            printf( " ] DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(0,i,N)]]] );
143            for( j=1; j<N; j++ )
144                printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,N)]]]  );  
145            printf( " ]\n" );  
146      }
147    }
148    
149    #endif
150    
151    #ifndef PASO_MPI
152    Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
153    #else
154    Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
155    #endif
156    {
157    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
158    index_t node0;    index_t node0;
159    Finley_Mesh* out;    Finley_Mesh* out;
# Line 41  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 165  Finley_Mesh* Finley_RectangularMesh_Hex8
165    N0=NE0+1;    N0=NE0+1;
166    N1=NE1+1;    N1=NE1+1;
167    N2=NE2+1;    N2=NE2+1;
168    #ifdef PASO_MPI
169      Paso_MPIInfo *mpi_info = NULL;
170    
171      /* get MPI information */
172      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
173      if (! Finley_noError())
174            return NULL;
175    #endif
176    
177    if (N0<=MIN(N1,N2)) {    if (N0<=MIN(N1,N2)) {
178       if (N1 <= N2) {       if (N1 <= N2) {
# Line 98  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 230  Finley_Mesh* Finley_RectangularMesh_Hex8
230    /*  allocate mesh: */    /*  allocate mesh: */
231        
232    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
233    /* TEMPFIX */  
234  #ifndef PASO_MPI  #ifndef PASO_MPI
235    out=Finley_Mesh_alloc(name,3,order);    out=Finley_Mesh_alloc(name,3,order);
236  #else  #else
237    out=Finley_Mesh_alloc(name,3,order,NULL);    out=Finley_Mesh_alloc(name,3,order,mpi_info);
238  #endif  #endif
239    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
240    
241  #ifndef PASO_MPI  #ifdef PASO_MPI
242      out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
243      if (useElementsOnFace) {
244         out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
245         out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
246      } else {
247         out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
248         out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
249      }
250      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
251    #else
252    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
253    if (useElementsOnFace) {    if (useElementsOnFace) {
254       out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);       out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
# Line 116  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 258  Finley_Mesh* Finley_RectangularMesh_Hex8
258       out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);       out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
259    }    }
260    out->Points=Finley_ElementFile_alloc(Point1,out->order);    out->Points=Finley_ElementFile_alloc(Point1,out->order);
 #else  
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,NULL);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,NULL);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,NULL);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,NULL);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,NULL);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,NULL);  
261  #endif  #endif
262    if (! Finley_noError()) {    if (! Finley_noError()) {
263        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
# Line 133  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 265  Finley_Mesh* Finley_RectangularMesh_Hex8
265    }    }
266    
267        
268    /*  allocate tables: */    /*  allocate tables: */
 #ifndef PASO_MPI    
269    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
270  #else  #ifdef PASO_MPI
271   /* TODO */    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
272  #endif  #endif
273    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
274    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
# Line 431  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 562  Finley_Mesh* Finley_RectangularMesh_Hex8
562        
563    /*  face elements done: */    /*  face elements done: */
564        
565    #ifdef PASO_MPI
566      /* make sure that the trivial distribution data is correct */
567        out->FaceElements->elementDistribution->numBoundary = 0;
568      out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = faceNECount;
569        out->Elements->elementDistribution->numBoundary = 0;
570      out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal = out->Elements->numElements;
571      out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
572      out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
573        
574        out->Nodes->degreeOfFreedomDistribution->numInternal = out->Nodes->degreeOfFreedomDistribution->numLocal;
575      out->Nodes->degreeOfFreedomDistribution->numBoundary = 0;
576    #endif
577    /*   condense the nodes: */    /*   condense the nodes: */
     
578    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
579    
580    #ifdef PASO_MPI
581      /* setup the CommBuffer */
582      Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
583      if ( !Finley_MPI_noError( mpi_info )) {
584        if( Finley_noError() )
585          Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
586        Paso_MPIInfo_dealloc( mpi_info );
587        Finley_Mesh_dealloc(out);
588        return NULL;
589      }
590    #endif
591    
592    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
593    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out) ;
594    
# Line 448  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 602  Finley_Mesh* Finley_RectangularMesh_Hex8
602    }    }
603    return out;    return out;
604  }  }
605    
606    #ifdef PASO_MPI
607    Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
608    {
609      dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
610      dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
611      bool_t dom_left, dom_right, dom_internal;
612    
613      index_t N0dom;
614      index_t firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1, node2;
615      index_t targetDomain=-1;
616      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
617      index_t *indexForward=NULL;
618      Finley_Mesh* out;
619    
620      char name[50];
621      Paso_MPIInfo *mpi_info = NULL;
622      double time0=Finley_timer();
623    
624      NE0=MAX(1,numElements[0]);
625      NE1=MAX(1,numElements[1]);
626      NE2=MAX(1,numElements[2]);
627      N0=NE0+1;
628      N1=NE1+1;
629      N2=NE2+1;
630    
631    
632      /* get MPI information */
633      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
634      if (! Finley_noError())
635            return NULL;
636    
637      /* use the serial version to generate the mesh for the 1-CPU case */
638      if( mpi_info->size==1 )
639      {
640        Paso_MPIInfo_dealloc( mpi_info );
641        out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace );
642        //print_mesh_statistics( out );
643            return out;
644      }    
645    
646      if( mpi_info->rank==0 )
647        domLeft = TRUE;
648      if( mpi_info->rank==mpi_info->size-1 )
649        domRight = TRUE;
650      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
651        domInternal = TRUE;
652    
653      /* dimensions of the local subdomain */
654      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );  
655    
656      /* count Degrees of Freedom along each axis */
657      NDOF0 = N0 - periodic[0];
658      NDOF1 = N1 - periodic[1];
659      NDOF2 = N2 - periodic[2];
660    
661      /* count face elements */
662      /* internal face elements */
663      NFaceElements = 0;
664      if( !periodic[0] )
665        NFaceElements += (domLeft+domRight)*NE1*NE2;
666      if( !periodic[1] )
667        NFaceElements += 2*numElementsInternal*NE2;
668      if( !periodic[2] )
669        NFaceElements += 2*numElementsInternal*NE1;
670      /* boundary face elements */
671        /* this is looks nasty, but it beats a bunch of nested if/then/else carry-on */
672      NFaceElements += 2*( 2 - (domLeft + domRight)*(!periodic[0]) )*( (!periodic[1])*NE2 + (!periodic[2])*NE1 );
673    
674      
675      /*  allocate mesh: */
676      sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
677    
678      out=Finley_Mesh_alloc(name,3,order,mpi_info);
679      if (! Finley_noError()) return NULL;
680    
681      out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
682      if (useElementsOnFace) {
683         out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
684         out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
685      } else {
686         out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
687         out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
688      }
689      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
690    
691      if (! Finley_noError()) {
692          Finley_Mesh_dealloc(out);
693          return NULL;
694      }
695    
696      
697      /*  allocate tables: */
698      Finley_NodeFile_allocTable(out->Nodes,(numNodesLocal+2-!periodic[0]*(domLeft+domRight))*N1*N2);
699      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, (DOFExternal[0]+DOFExternal[1])*NDOF1*NDOF2, 0 );
700      Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
701      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
702      if (! Finley_noError()) {
703          Finley_Mesh_dealloc(out);
704          return NULL;
705      }
706      
707      /*  set nodes: */
708      /* INTERNAL/BOUNDARY NODES */
709      #pragma omp parallel for private(i0,i1,i2,k)
710      for (k=0,i2=0;i2<N2;i2++) {
711        for (i1=0;i1<N1;i1++) {
712          for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {        
713            out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
714            out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
715            out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
716            out->Nodes->Id[k]=k;
717            out->Nodes->Tag[k]=0;
718            out->Nodes->degreeOfFreedom[k]=k;
719          }
720        }
721      }
722      if( domLeft && periodic[0] ) {
723        for (i2=0;i2<N2;i2++)  {
724          for (i1=0;i1<N1;i1++, k++) {
725            out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
726            out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
727            out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
728            out->Nodes->Id[k]=k;
729            out->Nodes->Tag[k]=0;
730            out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;
731          }
732        }
733        /* tag the faces for this special case */
734        if( !periodic[1] )
735        {
736          for( i2=0; i2<N2; i2++ ){
737            out->Nodes->Tag[k + (i2-N2)*N1     ] += 10;
738            out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;
739          }
740        }
741        if( !periodic[2] )
742        {
743          for( i1=0; i1<N1; i1++ ){
744            out->Nodes->Tag[k -N1*N2 +i1] += 100;
745            out->Nodes->Tag[k -N1    +i1] += 200;      
746          }
747        }
748      }
749      /* tags for the faces: */
750      N0dom = (numNodesLocal-periodicLocal[0]);
751      if (!periodic[2]) {    
752        for (i1=0;i1<N1;i1++) {
753          for (i0=0;i0<N0dom;i0++) {  
754             out->Nodes->Tag[i0 + N0dom*i1]+=100;
755             out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;
756           }
757         }
758      }
759      if (!periodic[1]) {
760        for (i2=0;i2<N2;i2++) {
761          for (i0=0;i0<N0dom;i0++) {
762             out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;
763             out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;
764          }
765        }
766      }
767      if (!periodic[0] && !domInternal ) {
768        for (i2=0;i2<N2;i2++) {
769          for (i1=0;i1<N1;i1++) {
770            if( domLeft )
771              out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;
772            if( domRight )
773              out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;
774          }
775        }
776      }
777        /* setup the forward communication data for the boundary nodes that we have just defined */
778        /* the case where there are 2 subdomains and periodic[0]=true has to be treated
779             as a special case to because the two domains have two interface boundaries to one-another */
780      indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );
781        if( mpi_info->size>2 || !periodic[0] ){
782            if( domInternal || domRight || periodicLocal[0] )
783            {
784                    for( int i=0; i<NDOF2; i++ )
785                        for( int j=0; j<NDOF1; j++ )
786                            indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
787                    targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
788                    Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
789            }
790            if( domInternal || domLeft || periodicLocal[1] )
791            {
792                for( int i=0; i<NDOF2; i++ )
793                    for( int j=0; j<NDOF1; j++ )
794                        indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
795                targetDomain = (mpi_info->rank+1) % mpi_info->size;
796                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
797            }
798        }
799        else {
800                for( int i=0; i<NDOF2; i++ )
801                    for( int j=0; j<NDOF1; j++ )
802                        indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;
803                targetDomain = (mpi_info->rank+1) % mpi_info->size;
804                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
805    
806                for( int i=0; i<NDOF2; i++ )
807                    for( int j=0; j<NDOF1; j++ )
808                        indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;
809                targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
810                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
811        }
812        
813      /* EXTERNAL NODES */
814      /* left hand boundary */
815      DOFcount = NDOF1*NDOF2*numDOFLocal;
816      if( (domLeft && periodic[0]) || !domLeft ) {
817        if( (domLeft && periodic[0]) )
818          for (i2=0;i2<N2;i2++)  {
819            for (i1=0;i1<N1;i1++, k++) {
820              out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];
821              out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
822              out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
823              out->Nodes->Id[k]=k;
824              out->Nodes->Tag[k]=0;
825                        out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
826            }
827          }
828        else
829          for (i2=0;i2<N2;i2++)  {
830            for (i1=0;i1<N1;i1++, k++) {
831              out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];
832              out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
833              out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
834              out->Nodes->Id[k]=k;
835              out->Nodes->Tag[k]=0;
836                        out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
837            }
838          }      
839            DOFcount += NDOF1*NDOF2;
840            targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
841            if( !periodic[1] ){
842                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
843            }
844            else {
845                for( int i=0; i<NDOF2; i++ )
846                    for( int j=0; j<NDOF1; j++ )
847                        indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
848                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
849            }
850                    
851        /* tag the faces for this special case */
852        if( !periodic[1] )
853        {
854          for( i1=0; i1<N1; i1++ ){
855            out->Nodes->Tag[k -N1*N2 +i1] += 10;
856            out->Nodes->Tag[k -N1    +i1] += 20;
857          }
858        }
859        if( periodic[2] )
860        {
861          for( i2=0; i2<N2; i2++ ){
862            out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;
863            out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;      
864          }
865        }
866      }
867      if( (domRight && periodic[0]) || !domRight )
868      {
869        if( domRight && periodic[0] )
870          for (i2=0;i2<N2;i2++)  {
871            for (i1=0;i1<N1;i1++, k++) {
872              out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];
873              out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
874              out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
875              out->Nodes->Id[k]=k;
876              out->Nodes->Tag[k]=0;
877                        out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
878            }
879          }
880        else
881          for (i2=0;i2<N2;i2++)  {
882            for (i1=0;i1<N1;i1++, k++) {
883              out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
884              out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
885              out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
886              out->Nodes->Id[k]=k;
887              out->Nodes->Tag[k]=0;
888                        out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;
889            }
890          }
891            DOFcount += NDOF1*NDOF2;
892    
893            targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;
894            if( !periodic[1] ){
895                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );
896            }
897            else {
898                for( int i=0; i<NDOF2; i++ )
899                    for( int j=0; j<NDOF1; j++ )
900                        indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;
901                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );
902            }
903                
904        /* tag the faces for this special case */
905        if( !periodic[1] )
906        {
907          for( i1=0; i1<N1; i1++ ){
908            out->Nodes->Tag[k -N1*N2 +i1] += 10;
909            out->Nodes->Tag[k -N1    +i1] += 20;
910          }
911        }
912        if( !periodic[2] )
913        {
914          for( i2=0; i2<N2; i2++ ){
915            out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;
916            out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;      
917          }
918        }
919      }
920        out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));
921        out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;
922      
923        TMPMEMFREE( indexForward );
924      /*   set the elements: */
925    
926      /* INTERNAL elements */
927      N0dom = (numNodesLocal-periodicLocal[0]);
928      k = 0;
929      #pragma omp parallel for private(i0,i1,i2,k,node0)
930      for (i2=0;i2<NE2;i2++) {
931        for (i1=0;i1<NE1;i1++) {
932          for (i0=0;i0<numElementsInternal;i0++,k++) {
933            node0=i0+i1*N0dom+N0dom*N1*i2;
934    
935            out->Elements->Id[k]=k;
936                    
937            out->Elements->Tag[k]=0;
938            out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
939    
940            out->Elements->Nodes[INDEX2(0,k,8)]=node0;
941            out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
942            out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;
943            out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;
944            out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;
945            out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;
946            out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;
947            out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;
948    
949          }
950        }
951      }
952        out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;
953        out->Elements->elementDistribution->numBoundary = 0;
954    
955      /* BOUNDARY Elements */
956      /* left boundary */
957      if( !domLeft )
958      {
959        for (i2=0;i2<NE2;i2++) {
960          node0 = numNodesLocal*N1*N2 + i2*N1;
961          for (i1=0;i1<NE1;i1++,node0++,k++) {
962            out->Elements->Id[k]=k;
963            out->Elements->Tag[k]=0;
964            out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
965    
966            out->Elements->Nodes[INDEX2(0,k,8)]=node0;
967            out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;
968            out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;
969            out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;
970            out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;
971            out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;
972            out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;
973            out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;    
974          }
975        }
976            out->Elements->elementDistribution->numBoundary += NE1*NE2;
977      }
978      /* the left periodic boundary is done a little differently to a left internal boundary */
979      else if( (domLeft && periodic[0]) )
980      {
981        for (i2=0;i2<NE2;i2++) {
982          node0 = numDOFLocal*N1*N2 + i2*N1;
983          node1 = node0 + N1*N2;
984          for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {
985    
986            out->Elements->Id[k]=k;
987            out->Elements->Tag[k]=0;
988            out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
989    
990            out->Elements->Nodes[INDEX2(0,k,8)]=node1;
991            out->Elements->Nodes[INDEX2(1,k,8)]=node0;
992            out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
993            out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;
994            out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;
995            out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
996            out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
997            out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;    
998          }
999        }  
1000            out->Elements->elementDistribution->numBoundary += NE1*NE2;
1001      }
1002      /* right boundary */
1003      if( !domRight || (domRight && periodic[0]) ){
1004        for (i2=0;i2<NE2;i2++) {
1005          for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {
1006                    node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;
1007                    node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;
1008    
1009            out->Elements->Id[k]=k;
1010            out->Elements->Tag[k]=0;
1011            out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
1012    
1013            out->Elements->Nodes[INDEX2(0,k,8)]=node1;
1014            out->Elements->Nodes[INDEX2(1,k,8)]=node0;
1015            out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;
1016            out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;
1017            out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;
1018            out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;
1019            out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;
1020            out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;  
1021          }
1022        }
1023            out->Elements->elementDistribution->numBoundary += NE1*NE2;
1024        }
1025    
1026        out->Elements->minColor=0;
1027      out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1028        out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;
1029        
1030      /*   face elements: */
1031      
1032      if  (useElementsOnFace) {
1033         NUMNODES=8;
1034      } else {
1035         NUMNODES=4;
1036      }
1037      totalNECount=k;
1038      faceNECount=0;
1039        idCount = totalNECount;
1040      
1041        /*   Do internal face elements for each boundary face first */
1042      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
1043        numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1044      if (!periodic[2]) {
1045         /*  elements on boundary 100 (x3=0): */
1046      
1047         #pragma omp parallel for private(i0,i1,k,node0)
1048         for (i1=0;i1<NE1;i1++) {
1049           for (i0=0; i0<numElementsInternal; i0++) {
1050             k=i0+numElementsInternal*i1+faceNECount;
1051             node0=i0+i1*numDOFLocal;
1052      
1053             out->FaceElements->Id[k]=idCount++;
1054             out->FaceElements->Tag[k]=100;
1055             out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);
1056    
1057             if  (useElementsOnFace) {
1058                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1059                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1060                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1061                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1062                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1063                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1064                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1065                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1066             } else {
1067                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1068                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1069                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;
1070                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1071             }
1072           }
1073         }
1074         totalNECount+=NE1*numElementsInternal;
1075         faceNECount+=NE1*numElementsInternal;
1076            
1077         /* **  elements on boundary 200 (x3=1) */
1078      
1079         #pragma omp parallel for private(i0,i1,k,node0)
1080         for (i1=0;i1<NE1;i1++) {
1081           for (i0=0;i0<numElementsInternal;i0++) {
1082             k=i0+numElementsInternal*i1+faceNECount;
1083             node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);
1084      
1085             out->FaceElements->Id[k]=idCount++;
1086             out->FaceElements->Tag[k]=200;
1087             out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;
1088    
1089             if  (useElementsOnFace) {
1090                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1091                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1092                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1093                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1094                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1095                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1096                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;
1097                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1098             } else {
1099                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;
1100                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;
1101                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;
1102                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;
1103             }
1104           }
1105         }
1106         totalNECount+=NE1*numElementsInternal;
1107         faceNECount+=NE1*numElementsInternal;
1108      }
1109      if (!periodic[0] && !domInternal) {
1110         /* **  elements on boundary 001 (x1=0): */
1111         if( domLeft ){
1112                 #pragma omp parallel for private(i1,i2,k,node0)
1113                 for (i2=0;i2<NE2;i2++) {
1114                     for (i1=0;i1<NE1;i1++) {
1115                         k=i1+NE1*i2+faceNECount;
1116                         node0=i1*numDOFLocal+numDOFLocal*N1*i2;
1117            
1118                         out->FaceElements->Id[k]=idCount++;
1119                         out->FaceElements->Tag[k]=1;
1120                         out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;
1121    
1122                         if  (useElementsOnFace) {
1123                                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1124                                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1125                                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1126                                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1127                                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;
1128                                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1129                                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1130                                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;
1131                         } else {
1132                                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1133                                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;
1134                                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1135                                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;
1136                         }
1137                     }
1138                 }
1139                 totalNECount+=NE1*NE2;
1140                 faceNECount+=NE1*NE2;
1141         }
1142         /* **  elements on boundary 002 (x1=1): */
1143             if( domRight ) {
1144                 #pragma omp parallel for private(i1,i2,k,node0)
1145                 for (i2=0;i2<NE2;i2++) {
1146                     for (i1=0;i1<NE1;i1++) {
1147                         k=i1+NE1*i2+faceNECount;
1148                         node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;
1149            
1150                         out->FaceElements->Id[k]=idCount++;
1151                         out->FaceElements->Tag[k]=2;
1152                         out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;
1153    
1154                         if  (useElementsOnFace) {
1155                                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1156                                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1157                                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1158                                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1159                                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1160                                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;
1161                                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1162                                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;
1163                         } else {
1164                                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
1165                                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;
1166                                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1167                                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1168                         }
1169                     }
1170                 }
1171           totalNECount+=NE1*NE2;
1172           faceNECount+=NE1*NE2;
1173         }
1174      }
1175      if (!periodic[1]) {
1176         /* **  elements on boundary 010 (x2=0): */
1177      
1178         #pragma omp parallel for private(i0,i2,k,node0)
1179         for (i2=0;i2<NE2;i2++) {
1180           for (i0=0;i0<numElementsInternal;i0++) {
1181             k=i0+numElementsInternal*i2+faceNECount;
1182             node0=i0+numDOFLocal*N1*i2;
1183      
1184             out->FaceElements->Id[k]=idCount++;
1185             out->FaceElements->Tag[k]=10;
1186             out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;
1187    
1188             if  (useElementsOnFace) {
1189                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1190                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1191                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1192                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1193                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1194                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;
1195                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;
1196                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1197             } else {
1198                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1199                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
1200                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;
1201                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1202             }
1203           }
1204         }
1205         totalNECount+=numElementsInternal*NE2;
1206         faceNECount+=numElementsInternal*NE2;
1207      
1208         /* **  elements on boundary 020 (x2=1): */
1209      
1210         #pragma omp parallel for private(i0,i2,k,node0)
1211         for (i2=0;i2<NE2;i2++) {
1212           for (i0=0;i0<numElementsInternal;i0++) {
1213             k=i0+numElementsInternal*i2+faceNECount;
1214             node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;
1215      
1216             out->FaceElements->Tag[k]=20;
1217             out->FaceElements->Id[k]=idCount++;
1218             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1219    
1220             if  (useElementsOnFace) {
1221                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1222                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1223                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1224                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1225                out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1226                out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;
1227                out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;
1228                out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;
1229             } else {
1230                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1231                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1232                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;
1233                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;
1234             }
1235    
1236           }
1237         }
1238         totalNECount+=numElementsInternal*NE2;
1239         faceNECount+=numElementsInternal*NE2;
1240      }
1241        out->FaceElements->elementDistribution->numInternal = faceNECount;
1242        
1243        /* now do the boundary face elements */
1244        /* LHS */
1245        if( !domLeft  )
1246        {
1247            if( !periodic[2] ) {
1248    
1249                /* x3=0 */
1250                for( i1=0; i1<NE1; i1++ )
1251                {
1252                    k = i1+faceNECount;
1253                    node0 = i1*numNodesLocal;
1254                    node1 = numNodesLocal*N1*N2 + i1;
1255    
1256                    out->FaceElements->Tag[k]=200;
1257                    out->FaceElements->Id[k]=idCount++;
1258                    out->FaceElements->Color[k]=0;
1259    
1260                    if( useElementsOnFace )
1261                    {
1262                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1263                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1264                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1265                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1266                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1267                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1268                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1269                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;
1270                    } else {
1271                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1272                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1273                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
1274                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1275                    }
1276                }  
1277                faceNECount  += NE1;
1278                totalNECount += NE1;
1279    
1280                /* x3=1 */
1281                for( i1=0; i1<NE1; i1++ )
1282                {
1283                    k = i1+faceNECount;
1284                    node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;
1285                    node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1286    
1287                    out->FaceElements->Tag[k]=200;
1288                    out->FaceElements->Id[k]=idCount++;
1289                    out->FaceElements->Color[k]=0;
1290    
1291                    if( useElementsOnFace )
1292                    {
1293                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1294                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1295                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1296                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1297                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1298                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1299                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;
1300                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1301                    } else {
1302                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1303                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;
1304                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1305                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1306                    }
1307                }  
1308                faceNECount  += NE1;
1309                totalNECount += NE1;
1310            }
1311    
1312            if( !periodic[1] ) {
1313                /* x2=0 */
1314                for( i2=0; i2<NE2; i2++ )
1315                {
1316                    k = i2+faceNECount;
1317                    node0 = i2*numNodesLocal*N1;
1318                    node1 = numNodesLocal*N1*N2 + i2*N1;
1319    
1320                    out->FaceElements->Tag[k]=20;
1321                    out->FaceElements->Id[k]=idCount++;
1322                    out->FaceElements->Color[k]=0;
1323    
1324                    if( useElementsOnFace )
1325                    {
1326                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1327                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1328                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1329                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1330                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1331                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;
1332                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1333                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1334                    } else {
1335                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1336                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1337                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;
1338                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1339                    }
1340                }  
1341                faceNECount  += NE2;
1342                totalNECount += NE2;
1343    
1344                /* x2=1 */
1345                for( i2=0; i2<NE2; i2++ )
1346                {
1347                    k = i2+faceNECount;
1348                    node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);
1349                    node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);
1350    
1351                    out->FaceElements->Tag[k]=20;
1352                    out->FaceElements->Id[k]=idCount++;
1353                    out->FaceElements->Color[k]=0;
1354    
1355                    if( useElementsOnFace )
1356                    {
1357                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1358                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1359                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;
1360                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1361                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1362                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1363                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1364                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1365                    } else {
1366                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1367                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
1368                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;
1369                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1370                    }
1371                }  
1372                faceNECount  += NE2;
1373                totalNECount += NE2;
1374            }
1375        }
1376      
1377        /* RHS */
1378        if( !domRight || periodicLocal[1] )
1379        {
1380            /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */
1381            if( domLeft && periodic[0] ){
1382                if( !periodic[2] ) {
1383    
1384                    /* x3=0 */
1385                    for( i1=0; i1<NE1; i1++ )
1386                    {
1387                        k = i1+faceNECount;
1388                        node0 = numDOFLocal*N1*N2 + i1;
1389                        node1 = numNodesLocal*N1*N2 + i1;
1390    
1391                        out->FaceElements->Tag[k]=200;
1392                        out->FaceElements->Id[k]=idCount++;
1393                        out->FaceElements->Color[k]=0;
1394    
1395                        if( useElementsOnFace )
1396                        {
1397                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1398                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1399                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1400                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1401                            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;
1402                            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;
1403                            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1404                            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;
1405                        } else {
1406                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1407                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;
1408                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
1409                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
1410                        }
1411                    }  
1412                    faceNECount  += NE1;
1413                    totalNECount += NE1;
1414    
1415                    /* x3=1 */
1416                    for( i1=0; i1<NE1; i1++ )
1417                    {
1418                        k = i1+faceNECount;
1419                        node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;
1420                        node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;
1421    
1422                        out->FaceElements->Tag[k]=200;
1423                        out->FaceElements->Id[k]=idCount++;
1424                        out->FaceElements->Color[k]=0;
1425                        
1426                        if( useElementsOnFace )
1427                        {
1428                            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1429                            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;
1430                            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;
1431                            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;
1432                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1433                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1434                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1435                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1436                        } else {
1437                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;
1438                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;
1439                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1440                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;
1441                        }
1442                    }  
1443                    faceNECount  += NE1;
1444                    totalNECount += NE1;
1445                }
1446    
1447                if( !periodic[1] ) {
1448                    /* x2=0 */
1449                    for( i2=0; i2<NE2; i2++ )
1450                    {
1451                        k = i2+faceNECount;
1452                        node0 = numDOFLocal*N1*N2 + i2*N1;
1453                        node1 = numNodesLocal*N1*N2 + i2*N1;
1454    
1455                        out->FaceElements->Tag[k]=20;
1456                        out->FaceElements->Id[k]=idCount++;
1457                        out->FaceElements->Color[k]=0;
1458    
1459                        if( useElementsOnFace )
1460                        {
1461                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1462                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1463                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1464                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1465                            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;
1466                            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;
1467                            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;
1468                            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;
1469                        } else {
1470                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;
1471                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
1472                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;
1473                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;
1474                        }
1475                    }  
1476                    faceNECount  += NE2;
1477                    totalNECount += NE2;
1478    
1479                    /* x2=1 */
1480                    for( i2=0; i2<NE2; i2++ )
1481                    {
1482                        k = i2+faceNECount;
1483                        node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;
1484                        node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;
1485    
1486                        out->FaceElements->Tag[k]=20;
1487                        out->FaceElements->Id[k]=idCount++;
1488                        out->FaceElements->Color[k]=0;
1489    
1490                        if( useElementsOnFace )
1491                        {
1492                            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;
1493                            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;
1494                            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;
1495                            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;
1496                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1497                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1498                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1499                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1500                        } else {
1501                            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;
1502                            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
1503                            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;
1504                            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;
1505                        }
1506                    }  
1507                    faceNECount  += NE2;
1508                    totalNECount += NE2;
1509                }
1510                
1511            }
1512            if( !periodic[2] ) {
1513                /* x3=0 */
1514                for( i1=0; i1<NE1; i1++ )
1515                {
1516                    k = i1+faceNECount;
1517                    node0 = numDOFLocal*(i1+1) - 1;
1518                    node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;
1519    
1520                    out->FaceElements->Tag[k]=200;
1521                    out->FaceElements->Id[k]=idCount++;
1522                    out->FaceElements->Color[k]=0;
1523    
1524                    if( useElementsOnFace )
1525                    {
1526                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1527                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1528                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1529                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1530                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;
1531                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1532                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1533                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;
1534                    } else {
1535                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1536                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;
1537                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;
1538                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;
1539                    }
1540                }  
1541                faceNECount  += NE1;
1542                totalNECount += NE1;
1543    
1544                /* x3=1 */
1545                for( i1=0; i1<NE1; i1++ )
1546                {
1547                    k = i1+faceNECount;
1548                    node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;
1549                    node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;
1550    
1551                    out->FaceElements->Tag[k]=200;
1552                    out->FaceElements->Id[k]=idCount++;
1553                    out->FaceElements->Color[k]=0;
1554                    
1555                    if( useElementsOnFace )
1556                    {
1557                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1558                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;
1559                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;
1560                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;
1561                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1562                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1563                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1564                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1565                    } else {
1566                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;
1567                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;
1568                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1569                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;
1570                    }
1571                }  
1572                faceNECount  += NE1;
1573                totalNECount += NE1;
1574            }
1575            if( !periodic[1] ) {
1576                /* x2=0 */
1577                for( i2=0; i2<NE2; i2++ )
1578                {
1579                    k = i2+faceNECount;
1580                    node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;
1581                    node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1582    
1583                    out->FaceElements->Tag[k]=20;
1584                    out->FaceElements->Id[k]=idCount++;
1585                    out->FaceElements->Color[k]=0;
1586    
1587                    if( useElementsOnFace )
1588                    {
1589                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1590                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1591                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1592                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1593                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;
1594                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;
1595                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;
1596                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1597                    } else {
1598                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
1599                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;
1600                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;
1601                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;
1602                    }
1603                }  
1604                faceNECount  += NE2;
1605                totalNECount += NE2;
1606    
1607                /* x2=1 */
1608                for( i2=0; i2<NE2; i2++ )
1609                {
1610                    k = i2+faceNECount;
1611                    node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;
1612                    node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;
1613    
1614                    out->FaceElements->Tag[k]=20;
1615                    out->FaceElements->Id[k]=idCount++;
1616                    out->FaceElements->Color[k]=0;
1617                    
1618                    if( useElementsOnFace ){
1619                        out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;
1620                        out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;
1621                        out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;
1622                        out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;
1623                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1624                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1625                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1626                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1627                    } else {
1628                        out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;
1629                        out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;
1630                        out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;
1631                        out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;
1632                    }
1633                }
1634                faceNECount  += NE2;
1635                totalNECount += NE2;
1636            }
1637        }
1638      out->FaceElements->minColor=0;
1639      out->FaceElements->maxColor=0;//23;
1640    
1641        out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;
1642      out->FaceElements->elementDistribution->numLocal = faceNECount;
1643    
1644    
1645      /* setup distribution info for other elements */
1646      out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;
1647      out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;
1648        
1649      /*   condense the nodes: */
1650      Finley_Mesh_resolveNodeIds( out );
1651    
1652      /* setup the CommBuffer */
1653      Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1654      if ( !Finley_MPI_noError( mpi_info )) {
1655        if( Finley_noError() )
1656          Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );
1657        Paso_MPIInfo_dealloc( mpi_info );
1658        Finley_Mesh_dealloc(out);
1659        return NULL;
1660      }
1661    
1662      Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );
1663    
1664      /* prepare mesh for further calculatuions:*/
1665      Finley_Mesh_prepare(out) ;
1666    
1667    //  print_mesh_statistics( out );
1668    
1669      #ifdef Finley_TRACE
1670      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1671      #endif
1672        
1673      if (! Finley_noError()) {
1674          Finley_Mesh_dealloc(out);
1675          return NULL;
1676      }
1677    
1678      return out;
1679    }
1680    #endif
1681    

Legend:
Removed from v.730  
changed lines
  Added in v.751

  ViewVC Help
Powered by ViewVC 1.1.26