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

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

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

revision 759 by bcumming, Thu Jun 29 01:53:23 2006 UTC revision 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC
# Line 112  static void domain_calculateDimension( i Line 112  static void domain_calculateDimension( i
112    DOFExternal[1] = nodesExternal[1];    DOFExternal[1] = nodesExternal[1];
113  }  }
114    
 void print_mesh_statistics( Finley_Mesh *out )  
 {  
   index_t i, j, N;  
     
   printf( "\nNodes\n=====\n\n" );  
     printf( "\t%d internal DOF\n\t%d boundary DOF\n\t%d local DOF\n\t%d external DOF\n", out->Nodes->degreeOfFreedomDistribution->numInternal, out->Nodes->degreeOfFreedomDistribution->numBoundary, out->Nodes->degreeOfFreedomDistribution->numLocal, out->Nodes->degreeOfFreedomDistribution->numExternal);  
   for( i=0; i<out->Nodes->numNodes; i++ )  
     printf( "node %d\t: id %d   \tDOF %d   \t: tag %d  \t: coordinates [%3g, %3g, %3g]\n", i, out->Nodes->Id[i], out->Nodes->degreeOfFreedom[i], out->Nodes->Tag[i], out->Nodes->Coordinates[INDEX2(0,i,3)], out->Nodes->Coordinates[INDEX2(1,i,3)], out->Nodes->Coordinates[INDEX2(2,i,3)] );  
   
   printf( "Elements\n========\n\n" );  
     printf( "\t%d internal\n\t%d boundary\n\t%d local\n", out->Elements->elementDistribution->numInternal, out->Elements->elementDistribution->numBoundary, out->Elements->elementDistribution->numLocal );  
     N = out->Elements->ReferenceElement->Type->numNodes;  
   for( i=0; i<out->Elements->numElements; i++ ){  
     printf( "element %d    \t: id %d  \t: nodes [ %3d, %3d, %3d, %3d, %3d, %3d, %3d, %3d ]", i, out->Elements->Id[i], out->Elements->Nodes[INDEX2(0,i,8)], out->Elements->Nodes[INDEX2(1,i,8)], out->Elements->Nodes[INDEX2(2,i,8)], out->Elements->Nodes[INDEX2(3,i,8)], out->Elements->Nodes[INDEX2(4,i,8)], out->Elements->Nodes[INDEX2(5,i,8)], out->Elements->Nodes[INDEX2(6,i,8)], out->Elements->Nodes[INDEX2(7,i,8)] );  
         printf( " DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(0,i,N)]]] );    
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->Elements->Nodes[INDEX2(j,i,N)]]]  );    
         printf( " ]\n" );    
   }  
   
     printf( "\nFace Elements\n==============\n\n" );  
     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 );  
     N = out->FaceElements->ReferenceElement->Type->numNodes;  
   for( i=0; i<out->FaceElements->numElements; i++ ){  
         printf( "face element %d \t: id %d  \t: nodes [ %3d", i, out->FaceElements->Id[i], out->FaceElements->Nodes[INDEX2(0,i,N)] );  
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->FaceElements->Nodes[INDEX2(j,i,N)]  );      
         printf( " ] DOF [ %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(0,i,N)]]] );  
         for( j=1; j<N; j++ )  
             printf( ", %3d", out->Nodes->degreeOfFreedom[out->Nodes->Id[out->FaceElements->Nodes[INDEX2(j,i,N)]]]  );    
         printf( " ]\n" );    
   }  
 }  
   
115  #endif  #endif
116    
117  #ifndef PASO_MPI  #ifdef PASO_MPI
118  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119  #else  #else
120  Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
121  #endif  #endif
122  {  {
123    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;
# Line 165  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 131  Finley_Mesh* Finley_RectangularMesh_Hex8
131    N0=NE0+1;    N0=NE0+1;
132    N1=NE1+1;    N1=NE1+1;
133    N2=NE2+1;    N2=NE2+1;
 #ifdef PASO_MPI  
   Paso_MPIInfo *mpi_info = NULL;  
   
   /* get MPI information */  
   mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  
   if (! Finley_noError())  
         return NULL;  
 #endif  
134    
135    if (N0<=MIN(N1,N2)) {    if (N0<=MIN(N1,N2)) {
136       if (N1 <= N2) {       if (N1 <= N2) {
# Line 290  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 248  Finley_Mesh* Finley_RectangularMesh_Hex8
248          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
249          out->Nodes->Tag[k]=0;          out->Nodes->Tag[k]=0;
250          out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);          out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
251    #ifdef PASO_MPI
252            out->Nodes->Dom[k]=NODE_INTERNAL;
253    #endif
254        }        }
255      }      }
256    }    }
# Line 332  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 293  Finley_Mesh* Finley_RectangularMesh_Hex8
293          out->Elements->Id[k]=k;          out->Elements->Id[k]=k;
294          out->Elements->Tag[k]=0;          out->Elements->Tag[k]=0;
295          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
296    #ifdef PASO_MPI
297            out->Elements->Dom[k]=ELEMENT_INTERNAL;
298    #endif
299    
300          out->Elements->Nodes[INDEX2(0,k,8)]=node0;          out->Elements->Nodes[INDEX2(0,k,8)]=node0;
301          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
# Line 371  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 335  Finley_Mesh* Finley_RectangularMesh_Hex8
335           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
336           out->FaceElements->Tag[k]=100;           out->FaceElements->Tag[k]=100;
337           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
338    #ifdef PASO_MPI
339            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
340    #endif
341    
342           if  (useElementsOnFace) {           if  (useElementsOnFace) {
343              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 403  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 370  Finley_Mesh* Finley_RectangularMesh_Hex8
370           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
371           out->FaceElements->Tag[k]=200;           out->FaceElements->Tag[k]=200;
372           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
373    #ifdef PASO_MPI
374            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
375    #endif
376    
377           if  (useElementsOnFace) {           if  (useElementsOnFace) {
378              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
# Line 438  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 408  Finley_Mesh* Finley_RectangularMesh_Hex8
408           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
409           out->FaceElements->Tag[k]=1;           out->FaceElements->Tag[k]=1;
410           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
411    #ifdef PASO_MPI
412            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
413    #endif
414    
415           if  (useElementsOnFace) {           if  (useElementsOnFace) {
416              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 470  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 443  Finley_Mesh* Finley_RectangularMesh_Hex8
443           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
444           out->FaceElements->Tag[k]=2;           out->FaceElements->Tag[k]=2;
445           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
446    #ifdef PASO_MPI
447            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
448    #endif
449    
450           if  (useElementsOnFace) {           if  (useElementsOnFace) {
451              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
# Line 503  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 479  Finley_Mesh* Finley_RectangularMesh_Hex8
479           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
480           out->FaceElements->Tag[k]=10;           out->FaceElements->Tag[k]=10;
481           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
482    #ifdef PASO_MPI
483            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
484    #endif
485    
486           if  (useElementsOnFace) {           if  (useElementsOnFace) {
487              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 535  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 514  Finley_Mesh* Finley_RectangularMesh_Hex8
514           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
515           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
516           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
517    #ifdef PASO_MPI
518            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
519    #endif
520    
521           if  (useElementsOnFace) {           if  (useElementsOnFace) {
522              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
# Line 560  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 542  Finley_Mesh* Finley_RectangularMesh_Hex8
542    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
543    out->FaceElements->maxColor=23;    out->FaceElements->maxColor=23;
544        
   /*  face elements done: */  
     
545  #ifdef PASO_MPI  #ifdef PASO_MPI
546    /* make sure that the trivial distribution data is correct */      Finley_ElementFile_setDomainFlags( out->Elements );
547      out->FaceElements->elementDistribution->numBoundary = 0;      Finley_ElementFile_setDomainFlags( out->FaceElements );
548    out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal = faceNECount;      Finley_ElementFile_setDomainFlags( out->ContactElements );
549      out->Elements->elementDistribution->numBoundary = 0;      Finley_ElementFile_setDomainFlags( out->Points );
550    out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal = out->Elements->numElements;  
551    out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;      /* reorder the degrees of freedom */
552    out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;      Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
       
     out->Nodes->degreeOfFreedomDistribution->numInternal = out->Nodes->degreeOfFreedomDistribution->numLocal;  
   out->Nodes->degreeOfFreedomDistribution->numBoundary = 0;  
553  #endif  #endif
554        
555    /*   condense the nodes: */    /*   condense the nodes: */
556    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
557    
 #ifdef PASO_MPI  
   /* setup the CommBuffer */  
   Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
   if ( !Finley_MPI_noError( mpi_info )) {  
     if( Finley_noError() )  
       Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
 #endif  
   
558    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
559    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out) ;
560    
# Line 606  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 572  Finley_Mesh* Finley_RectangularMesh_Hex8
572  #ifdef PASO_MPI  #ifdef PASO_MPI
573  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
574  {  {
575    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,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
576    dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;    dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
577    bool_t dom_left, dom_right, dom_internal;    bool_t dom_left, dom_right, dom_internal;
578      index_t firstNode=0, DOFcount=0, node0, node1, node2;
579    index_t N0dom;    index_t targetDomain=-1, firstNodeConstruct, j;
580    index_t firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1, node2;    bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
581    index_t targetDomain=-1;      index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;  
   index_t *indexForward=NULL;  
582    Finley_Mesh* out;    Finley_Mesh* out;
583    
584    char name[50];    char name[50];
585    Paso_MPIInfo *mpi_info = NULL;    Paso_MPIInfo *mpi_info = NULL;
586    double time0=Finley_timer();    double time0=Finley_timer();
587    
588        index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
589        index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
590        index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
591        index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
592        index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
593        index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
594    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
595    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
596    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
# Line 637  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 607  Finley_Mesh* Finley_RectangularMesh_Hex8
607    /* use the serial version to generate the mesh for the 1-CPU case */    /* use the serial version to generate the mesh for the 1-CPU case */
608    if( mpi_info->size==1 )    if( mpi_info->size==1 )
609    {    {
610      Paso_MPIInfo_dealloc( mpi_info );      out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info );
     out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace );  
     //print_mesh_statistics( out );  
611          return out;          return out;
612    }        }    
613    
# Line 664  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 632  Finley_Mesh* Finley_RectangularMesh_Hex8
632    if( !periodic[0] )    if( !periodic[0] )
633      NFaceElements += (domLeft+domRight)*NE1*NE2;      NFaceElements += (domLeft+domRight)*NE1*NE2;
634    if( !periodic[1] )    if( !periodic[1] )
635      NFaceElements += 2*numElementsInternal*NE2;      NFaceElements += 2*numElementsLocal*NE2;
636    if( !periodic[2] )    if( !periodic[2] )
637      NFaceElements += 2*numElementsInternal*NE1;      NFaceElements += 2*numElementsLocal*NE1;
   /* boundary face elements */  
     /* this is looks nasty, but it beats a bunch of nested if/then/else carry-on */  
   NFaceElements += 2*( 2 - (domLeft + domRight)*(!periodic[0]) )*( (!periodic[1])*NE2 + (!periodic[2])*NE1 );  
   
638        
639        boundaryLeft = !domLeft || periodicLocal[0];
640        boundaryRight = !domRight || periodicLocal[1];
641        N0t = numNodesLocal + boundaryRight + boundaryLeft;
642        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
643        firstNodeConstruct = firstNode - boundaryLeft;
644        firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
645    
646    /*  allocate mesh: */    /*  allocate mesh: */
647    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
648    
# Line 692  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 663  Finley_Mesh* Finley_RectangularMesh_Hex8
663        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
664        return NULL;        return NULL;
665    }    }
   
666        
667    /*  allocate tables: */    /*  allocate tables: */
668    Finley_NodeFile_allocTable(out->Nodes,(numNodesLocal+2-!periodic[0]*(domLeft+domRight))*N1*N2);    Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
669    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, (DOFExternal[0]+DOFExternal[1])*NDOF1*NDOF2, 0 );    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
670    Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);    Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
671    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
672    if (! Finley_noError()) {    if (! Finley_noError()) {
# Line 710  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 680  Finley_Mesh* Finley_RectangularMesh_Hex8
680    #pragma omp parallel for private(i0,i1,i2,k)    #pragma omp parallel for private(i0,i1,i2,k)
681    for (i2=0;i2<N2;i2++) {    for (i2=0;i2<N2;i2++) {
682      for (i1=0;i1<N1;i1++) {      for (i1=0;i1<N1;i1++) {
683        for (i0=0;i0<numNodesLocal-domLeft*periodic[0];i0++,k++) {                for (i0=0;i0<N0t;i0++,k++) {        
684          out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];          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]=k;  
       }  
     }  
   }  
   if( domLeft && periodic[0] ) {  
     for (i2=0;i2<N2;i2++)  {  
       for (i1=0;i1<N1;i1++, k++) {  
         out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];  
685          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
686          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
687          out->Nodes->Id[k]=k;          out->Nodes->Id[k]=k;
688          out->Nodes->Tag[k]=0;          out->Nodes->Tag[k]=0;
689          out->Nodes->degreeOfFreedom[k]=i1*numNodesLocal + i2*numNodesLocal*N1;          out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
690        }                  out->Nodes->Dom[k]=NODE_INTERNAL;
     }  
     /* tag the faces for this special case */  
     if( !periodic[1] )  
     {  
       for( i2=0; i2<N2; i2++ ){  
         out->Nodes->Tag[k + (i2-N2)*N1     ] += 10;  
         out->Nodes->Tag[k + (i2+1-N2)*N1 -1] += 20;  
       }  
     }  
     if( !periodic[2] )  
     {  
       for( i1=0; i1<N1; i1++ ){  
         out->Nodes->Tag[k -N1*N2 +i1] += 100;  
         out->Nodes->Tag[k -N1    +i1] += 200;        
691        }        }
692      }      }
693    }    }
694    /* tags for the faces: */  
695    N0dom = (numNodesLocal-periodicLocal[0]);      /* mark the nodes that reference external and boundary DOF as such */
696        if( boundaryLeft ){
697            for( i1=0; i1<N1; i1++ )
698                for( i2=0; i2<N2; i2++ ) {
699                    out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
700                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
701                }
702        }
703        if( boundaryRight ){
704            for( i1=0; i1<N1; i1++ )
705                for( i2=0; i2<N2; i2++ ) {
706                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
707                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
708                }
709        }
710        if( periodicLocal[0] ){
711            for( i1=0; i1<N1; i1++ )
712                for( i2=0; i2<N2; i2++ ) {
713                    out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
714                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
715                }
716        }
717            
718      /* tag Nodes that are referenced by face elements */
719    if (!periodic[2]) {        if (!periodic[2]) {    
720      for (i1=0;i1<N1;i1++) {      for (i1=0;i1<N1;i1++) {
721        for (i0=0;i0<N0dom;i0++) {          for (i0=0;i0<N0t;i0++) {  
722           out->Nodes->Tag[i0 + N0dom*i1]+=100;           out->Nodes->Tag[i0 + N0t*i1]+=100;
723           out->Nodes->Tag[i0 + N0dom*i1 + N0dom*N1*(N2-1)]+=200;           out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
724         }         }
725       }       }
726    }    }
727    if (!periodic[1]) {    if (!periodic[1]) {
728      for (i2=0;i2<N2;i2++) {      for (i2=0;i2<N2;i2++) {
729        for (i0=0;i0<N0dom;i0++) {        for (i0=0;i0<N0t;i0++) {
730           out->Nodes->Tag[i0 + i2*N1*N0dom]+=10;           out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
731           out->Nodes->Tag[i0 + (i2+1)*N1*N0dom-N0dom]+=20;           out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
732        }        }
733      }      }
734    }    }
# Line 769  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 736  Finley_Mesh* Finley_RectangularMesh_Hex8
736      for (i2=0;i2<N2;i2++) {      for (i2=0;i2<N2;i2++) {
737        for (i1=0;i1<N1;i1++) {        for (i1=0;i1<N1;i1++) {
738          if( domLeft )          if( domLeft )
739            out->Nodes->Tag[i1*N0dom + i2*N0dom*N1]+=1;            out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
740          if( domRight )          if( domRight )
741            out->Nodes->Tag[(i1+1)*N0dom-1 + i2*N0dom*N1]+=2;            out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
742        }        }
743      }      }
744    }    }
     /* setup the forward communication data for the boundary nodes that we have just defined */  
     /* the case where there are 2 subdomains and periodic[0]=true has to be treated  
          as a special case to because the two domains have two interface boundaries to one-another */  
   indexForward = TMPMEMALLOC( NDOF1*NDOF2, index_t );  
     if( mpi_info->size>2 || !periodic[0] ){  
         if( domInternal || domRight || periodicLocal[0] )  
         {  
                 for( int i=0; i<NDOF2; i++ )  
                     for( int j=0; j<NDOF1; j++ )  
                         indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;  
                 targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;  
                 Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
         }  
         if( domInternal || domLeft || periodicLocal[1] )  
         {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;  
             targetDomain = (mpi_info->rank+1) % mpi_info->size;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
         }  
     }  
     else {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*(j+1)-1+NDOF1*numDOFLocal*i;  
             targetDomain = (mpi_info->rank+1) % mpi_info->size;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
   
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = numDOFLocal*j+NDOF1*numDOFLocal*i;  
             targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;  
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
     }  
       
   /* EXTERNAL NODES */  
   /* left hand boundary */  
   DOFcount = NDOF1*NDOF2*numDOFLocal;  
   if( (domLeft && periodic[0]) || !domLeft ) {  
     if( (domLeft && periodic[0]) )  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=(1.-DBLE(1)/DBLE(N0-1))*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }  
     else  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=(DBLE(firstNode-1)/DBLE(N0-1))*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }        
         DOFcount += NDOF1*NDOF2;  
         targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;  
         if( !periodic[1] ){  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );  
         }  
         else {  
             for( int i=0; i<NDOF2; i++ )  
                 for( int j=0; j<NDOF1; j++ )  
                     indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;  
             Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );  
         }  
                   
     /* tag the faces for this special case */  
     if( !periodic[1] )  
     {  
       for( i1=0; i1<N1; i1++ ){  
         out->Nodes->Tag[k -N1*N2 +i1] += 10;  
         out->Nodes->Tag[k -N1    +i1] += 20;  
       }  
     }  
     if( periodic[2] )  
     {  
       for( i2=0; i2<N2; i2++ ){  
         out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;  
         out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;        
       }  
     }  
   }  
   if( (domRight && periodic[0]) || !domRight )  
   {  
     if( domRight && periodic[0] )  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }  
     else  
       for (i2=0;i2<N2;i2++)  {  
         for (i1=0;i1<N1;i1++, k++) {  
           out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];  
           out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
           out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];  
           out->Nodes->Id[k]=k;  
           out->Nodes->Tag[k]=0;  
                     out->Nodes->degreeOfFreedom[k]=DOFcount+i1%NDOF1+(i2%NDOF2)*NDOF1;  
         }  
       }  
         DOFcount += NDOF1*NDOF2;  
745    
746          targetDomain = mpi_info->rank+1 < mpi_info->size? mpi_info->rank+1 : 0;      /* form the boudary communication information */
747          if( !periodic[1] ){      forwardDOF  = MEMALLOC(NDOF1*NDOF2,index_t);
748              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, out->Nodes->degreeOfFreedom + (k-NDOF1*NDOF2) );      backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
749          }      if( !(mpi_info->size==2 && periodicLocal[0])){
750          else {          if( boundaryLeft  ) {
751              for( int i=0; i<NDOF2; i++ )              targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
752                  for( int j=0; j<NDOF1; j++ )              for( i2=0; i2<NDOF2; i2++ ){
753                      indexForward[j+i*NDOF1] = DOFcount - NDOF1*NDOF2 + j + i*NDOF1;                  for( i1=0; i1<NDOF1; i1++ ){
754              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, indexForward );                      forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
755          }                      backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
756                                }
757      /* tag the faces for this special case */              }
758      if( !periodic[1] )              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
759      {              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
760        for( i1=0; i1<N1; i1++ ){          }
761          out->Nodes->Tag[k -N1*N2 +i1] += 10;          if( boundaryRight ) {
762          out->Nodes->Tag[k -N1    +i1] += 20;              targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
763        }              for( i2=0; i2<NDOF2; i2++ ){
764      }                  for( i1=0; i1<NDOF1; i1++ ){
765      if( !periodic[2] )                      forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
766      {                      backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
767        for( i2=0; i2<N2; i2++ ){                  }
768          out->Nodes->Tag[k +(i2-N2)*N1     ] += 100;              }
769          out->Nodes->Tag[k +(i2-N2+1)*N1 -1] += 200;                    Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
770        }              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
771      }          }
772    }      } else{
773      out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*NDOF2*(numDOFLocal - 2 + domRight*(!periodic[0]) + domLeft*(!periodic[0]));          /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
774      out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal - out->Nodes->degreeOfFreedomDistribution->numInternal;          targetDomain = 1;
775              
776      TMPMEMFREE( indexForward );          for( i2=0; i2<NDOF2; i2++ ){
777                for( i1=0; i1<NDOF1; i1++ ){
778                    forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
779                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
780                }
781            }
782            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
783            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
784    
785            for( i2=0; i2<NDOF2; i2++ ){
786                for( i1=0; i1<NDOF1; i1++ ){
787                    forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
788                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
789                }
790            }
791            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
792            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
793        }
794        MEMFREE( forwardDOF );
795        MEMFREE( backwardDOF );
796    /*   set the elements: */    /*   set the elements: */
797    
798    /* INTERNAL elements */    /* INTERNAL elements */
   N0dom = (numNodesLocal-periodicLocal[0]);  
799    k = 0;    k = 0;
800    #pragma omp parallel for private(i0,i1,i2,k,node0)    #pragma omp parallel for private(i0,i1,i2,k,node0)
801    for (i2=0;i2<NE2;i2++) {    for (i2=0;i2<NE2;i2++) {
802      for (i1=0;i1<NE1;i1++) {      for (i1=0;i1<NE1;i1++) {
803        for (i0=0;i0<numElementsInternal;i0++,k++) {        for (i0=0;i0<numElementsLocal;i0++,k++) {
804          node0=i0+i1*N0dom+N0dom*N1*i2;                  node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t :  i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
   
805          out->Elements->Id[k]=k;          out->Elements->Id[k]=k;
806                                    
807          out->Elements->Tag[k]=0;          out->Elements->Tag[k]=0;
808          out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
809                    out->Elements->Dom[k]=ELEMENT_INTERNAL;
810    
811          out->Elements->Nodes[INDEX2(0,k,8)]=node0;          out->Elements->Nodes[INDEX2(0,k,8)]=node0;
812          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
813          out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0dom+1;          out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
814          out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0dom;          out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
815          out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0dom*N1;          out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
816          out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0dom*N1+1;          out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
817          out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0dom*N1+N0dom+1;          out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
818          out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0dom*N1+N0dom;          out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
819    
820        }        }
821      }      }
822    }    }
     out->Elements->elementDistribution->numInternal = NE1*NE2*numElementsInternal;  
     out->Elements->elementDistribution->numBoundary = 0;  
   
   /* BOUNDARY Elements */  
   /* left boundary */  
   if( !domLeft )  
   {  
     for (i2=0;i2<NE2;i2++) {  
       node0 = numNodesLocal*N1*N2 + i2*N1;  
       for (i1=0;i1<NE1;i1++,node0++,k++) {  
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,8)]=i1*numNodesLocal + i2*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(2,k,8)]=(i1+1)*numNodesLocal + i2*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=i1*numNodesLocal + (i2+1)*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=(i1+1)*numNodesLocal + (i2+1)*numNodesLocal*N1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N1+1;      
       }  
     }  
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
   }  
   /* the left periodic boundary is done a little differently to a left internal boundary */  
   else if( (domLeft && periodic[0]) )  
   {  
     for (i2=0;i2<NE2;i2++) {  
       node0 = numDOFLocal*N1*N2 + i2*N1;  
       node1 = node0 + N1*N2;  
       for (i1=0;i1<NE1;i1++,k++,node0++,node1++) {  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node1+1;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node1+N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node1+N1+1;      
       }  
     }    
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
   }  
   /* right boundary */  
   if( !domRight || (domRight && periodic[0]) ){  
     for (i2=0;i2<NE2;i2++) {  
       for (i1=0;i1<NE1;i1++,node0++,node1+=numDOFLocal,k++) {  
                 node1 = numDOFLocal -1 + numDOFLocal*i1 + N1*numDOFLocal*i2;  
                 node0 = (numNodesLocal + domInternal + periodicLocal[0])*N1*N2 + i2*N1 + i1;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0;//COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node1+N0dom;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node1+N0dom*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N1+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node1+N0dom*N1+N0dom;    
       }  
     }  
         out->Elements->elementDistribution->numBoundary += NE1*NE2;  
     }  
   
823      out->Elements->minColor=0;      out->Elements->minColor=0;
824    out->Elements->maxColor=0;//COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);    out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
825      out->Elements->elementDistribution->numLocal = out->Elements->elementDistribution->numInternal + out->Elements->elementDistribution->numBoundary;      if( boundaryLeft )
826            for( i2=0; i2<NE2; i2++ )
827                for( i1=0; i1<NE1; i1++ )
828                    out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
829        if( boundaryRight )
830            for( i2=0; i2<NE2; i2++ )
831                for( i1=0; i1<NE1; i1++ )
832                    out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
833            
834        Finley_ElementFile_setDomainFlags( out->Elements );
835            
836    /*   face elements: */    /*   face elements: */
     
837    if  (useElementsOnFace) {    if  (useElementsOnFace) {
838       NUMNODES=8;       NUMNODES=8;
839    } else {    } else {
840       NUMNODES=4;       NUMNODES=4;
841    }    }
842    totalNECount=k;    totalNECount=out->Elements->numElements;
843    faceNECount=0;    faceNECount=0;
844      idCount = totalNECount;      idCount = totalNECount;
845        
     /*   Do internal face elements for each boundary face first */  
846    /*   these are the quadrilateral elements on boundary 1 (x3=0): */    /*   these are the quadrilateral elements on boundary 1 (x3=0): */
847      numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];      numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
848    if (!periodic[2]) {    if (!periodic[2]) {
849       /*  elements on boundary 100 (x3=0): */       /*  elements on boundary 100 (x3=0): */
850        
851       #pragma omp parallel for private(i0,i1,k,node0)       #pragma omp parallel for private(i0,i1,k)
852       for (i1=0;i1<NE1;i1++) {       for (i1=0;i1<NE1;i1++) {
853         for (i0=0; i0<numElementsInternal; i0++) {         for (i0=0; i0<numElementsLocal; i0++) {
854           k=i0+numElementsInternal*i1+faceNECount;           k=i0+numElementsLocal*i1+faceNECount;
855           node0=i0+i1*numDOFLocal;                   kk=i0 + i1*numElementsLocal;
856                     facePerm = face2;
857        
858           out->FaceElements->Id[k]=idCount++;           out->FaceElements->Id[k]=idCount++;
859           out->FaceElements->Tag[k]=100;           out->FaceElements->Tag[k]=100;
860           out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2);           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
861             out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
862    
863           if  (useElementsOnFace) {                   for( j=0; j<NUMNODES; j++ )
864              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
          }  
865         }         }
866       }       }
867       totalNECount+=NE1*numElementsInternal;         if( boundaryLeft ){
868       faceNECount+=NE1*numElementsInternal;              for( i1=0; i1<NE1; i1++ )
869                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
870                if( periodicLocal[0] )
871                    for( i1=0; i1<NE1; i1++ )
872                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
873             }
874           if( boundaryRight )
875                for( i1=0; i1<NE1; i1++ )
876                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
877         totalNECount+=NE1*numElementsLocal;
878         faceNECount+=NE1*numElementsLocal;
879                    
880       /* **  elements on boundary 200 (x3=1) */       /* **  elements on boundary 200 (x3=1) */
881        
882       #pragma omp parallel for private(i0,i1,k,node0)       #pragma omp parallel for private(i0,i1,k)
883       for (i1=0;i1<NE1;i1++) {       for (i1=0;i1<NE1;i1++) {
884         for (i0=0;i0<numElementsInternal;i0++) {         for (i0=0;i0<numElementsLocal;i0++) {
885           k=i0+numElementsInternal*i1+faceNECount;           k=i0+numElementsLocal*i1+faceNECount;
886           node0=i0+i1*numDOFLocal+numDOFLocal*N1*(NE2-1);                   kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
887                     facePerm = face3;
888        
889           out->FaceElements->Id[k]=idCount++;           out->FaceElements->Id[k]=idCount++;
890           out->FaceElements->Tag[k]=200;           out->FaceElements->Tag[k]=200;
891           out->FaceElements->Color[k]=0;//(i0%2)+2*(i1%2)+4;           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
892             out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
893    
894           if  (useElementsOnFace) {                   for( j=0; j<NUMNODES; j++ )
895              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ numDOFLocal * N1;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ numDOFLocal * N1+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ numDOFLocal * N1+numDOFLocal;  
          }  
896         }         }
897       }       }
898       totalNECount+=NE1*numElementsInternal;         if( boundaryLeft ){
899       faceNECount+=NE1*numElementsInternal;              for( i1=0; i1<NE1; i1++ )
900                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
901                if( periodicLocal[0] )
902                    for( i1=0; i1<NE1; i1++ )
903                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
904             }
905           if( boundaryRight )
906                for( i1=0; i1<NE1; i1++ )
907                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
908         totalNECount+=NE1*numElementsLocal;
909         faceNECount+=NE1*numElementsLocal;
910    }    }
911    if (!periodic[0] && !domInternal) {    if (!periodic[0] && !domInternal) {
912       /* **  elements on boundary 001 (x1=0): */       /* **  elements on boundary 001 (x1=0): */
913       if( domLeft ){       if( domLeft ){
914               #pragma omp parallel for private(i1,i2,k,node0)               #pragma omp parallel for private(i1,i2,k)
915               for (i2=0;i2<NE2;i2++) {               for (i2=0;i2<NE2;i2++) {
916                   for (i1=0;i1<NE1;i1++) {                   for (i1=0;i1<NE1;i1++) {
917                       k=i1+NE1*i2+faceNECount;                       k=i1+NE1*i2+faceNECount;
918                       node0=i1*numDOFLocal+numDOFLocal*N1*i2;                       kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
919                         facePerm = face0;
920                    
921                       out->FaceElements->Id[k]=idCount++;                       out->FaceElements->Id[k]=idCount++;
922                       out->FaceElements->Tag[k]=1;                       out->FaceElements->Tag[k]=1;
923                       out->FaceElements->Color[k]=0; //(i2%2)+2*(i1%2)+8;                       out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
924                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
925    
926                       if  (useElementsOnFace) {                       for( j=0; j<NUMNODES; j++ )
927                              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                          out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
                             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal+1;  
                      } else {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal;  
                      }  
928                   }                   }
929               }               }
930               totalNECount+=NE1*NE2;               totalNECount+=NE1*NE2;
# Line 1142  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 932  Finley_Mesh* Finley_RectangularMesh_Hex8
932       }       }
933       /* **  elements on boundary 002 (x1=1): */       /* **  elements on boundary 002 (x1=1): */
934           if( domRight ) {           if( domRight ) {
935               #pragma omp parallel for private(i1,i2,k,node0)               #pragma omp parallel for private(i1,i2,k)
936               for (i2=0;i2<NE2;i2++) {               for (i2=0;i2<NE2;i2++) {
937                   for (i1=0;i1<NE1;i1++) {                   for (i1=0;i1<NE1;i1++) {
938                       k=i1+NE1*i2+faceNECount;                       k=i1+NE1*i2+faceNECount;
939                       node0=(numDOFLocal-2)+i1*numDOFLocal+numDOFLocal*N1*i2 ;                       kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
940                         facePerm = face1;
941                    
942                       out->FaceElements->Id[k]=idCount++;                       out->FaceElements->Id[k]=idCount++;
943                       out->FaceElements->Tag[k]=2;                       out->FaceElements->Tag[k]=2;
944                       out->FaceElements->Color[k]=0;//(i2%2)+2*(i1%2)+12;               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
945                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
946    
947                       if  (useElementsOnFace) {                       for( j=0; j<NUMNODES; j++ )
948                              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                          out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal*N1;  
                      } else {  
                             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
                             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
                      }  
949                   }                   }
950               }               }
951         totalNECount+=NE1*NE2;         totalNECount+=NE1*NE2;
# Line 1176  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 955  Finley_Mesh* Finley_RectangularMesh_Hex8
955    if (!periodic[1]) {    if (!periodic[1]) {
956       /* **  elements on boundary 010 (x2=0): */       /* **  elements on boundary 010 (x2=0): */
957        
958       #pragma omp parallel for private(i0,i2,k,node0)       #pragma omp parallel for private(i0,i2,k)
959       for (i2=0;i2<NE2;i2++) {       for (i2=0;i2<NE2;i2++) {
960         for (i0=0;i0<numElementsInternal;i0++) {         for (i0=0;i0<numElementsLocal;i0++) {
961           k=i0+numElementsInternal*i2+faceNECount;           k=i0+numElementsLocal*i2+faceNECount;
962           node0=i0+numDOFLocal*N1*i2;                   kk=i0+numElementsLocal*NE1*i2;
963                     facePerm = face4;
964        
965           out->FaceElements->Id[k]=idCount++;           out->FaceElements->Id[k]=idCount++;
966           out->FaceElements->Tag[k]=10;           out->FaceElements->Tag[k]=10;
967           out->FaceElements->Color[k]=0;//(i0%2)+2*(i2%2)+16;           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
968             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
969    
970           if  (useElementsOnFace) {                   for( j=0; j<NUMNODES; j++ )
971              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
          }  
972         }         }
973       }       }
974       totalNECount+=numElementsInternal*NE2;         if( boundaryLeft ){
975       faceNECount+=numElementsInternal*NE2;              for( i2=0; i2<NE2; i2++ )
976                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
977                if( periodicLocal[0] )
978                    for( i2=0; i2<NE2; i2++ )
979                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
980             }
981           if( boundaryRight )
982                for( i2=0; i2<NE2; i2++ )
983                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
984         totalNECount+=numElementsLocal*NE2;
985         faceNECount+=numElementsLocal*NE2;
986        
987       /* **  elements on boundary 020 (x2=1): */       /* **  elements on boundary 020 (x2=1): */
988        
989       #pragma omp parallel for private(i0,i2,k,node0)       #pragma omp parallel for private(i0,i2,k)
990       for (i2=0;i2<NE2;i2++) {       for (i2=0;i2<NE2;i2++) {
991         for (i0=0;i0<numElementsInternal;i0++) {         for (i0=0;i0<numElementsLocal;i0++) {
992           k=i0+numElementsInternal*i2+faceNECount;           k=i0+numElementsLocal*i2+faceNECount;
993           node0=i0+(NE1-1)*numDOFLocal+numDOFLocal*N1*i2;                   kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
994                     facePerm = face5;
995        
996           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
997           out->FaceElements->Id[k]=idCount++;           out->FaceElements->Id[k]=idCount++;
998             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
999           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1000    
1001           if  (useElementsOnFace) {                   for( j=0; j<NUMNODES; j++ )
1002              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numDOFLocal*N1+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal+1;  
          }  
   
1003         }         }
1004       }       }
1005       totalNECount+=numElementsInternal*NE2;         if( boundaryLeft ){
1006       faceNECount+=numElementsInternal*NE2;              for( i2=0; i2<NE2; i2++ )
1007                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1008                if( periodicLocal[0] )
1009                    for( i2=0; i2<NE2; i2++ )
1010                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1011             }
1012           if( boundaryRight )
1013                for( i2=0; i2<NE2; i2++ )
1014                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1015         totalNECount+=numElementsLocal*NE2;
1016         faceNECount+=numElementsLocal*NE2;
1017    }    }
1018      out->FaceElements->elementDistribution->numInternal = faceNECount;      out->FaceElements->elementDistribution->numInternal = faceNECount;
1019            
     /* now do the boundary face elements */  
     /* LHS */  
     if( !domLeft  )  
     {  
         if( !periodic[2] ) {  
   
             /* x3=0 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = i1*numNodesLocal;  
                 node1 = numNodesLocal*N1*N2 + i1;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numNodesLocal*N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
   
             /* x3=1 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numNodesLocal*N1*(NE2-1) + i1*numNodesLocal;  
                 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
         }  
   
         if( !periodic[1] ) {  
             /* x2=0 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = i2*numNodesLocal*N1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
   
             /* x2=1 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = i2*numNodesLocal*N1 + numNodesLocal*(NE1-1);  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1);  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+numNodesLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal*N1+numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
         }  
     }  
     
     /* RHS */  
     if( !domRight || periodicLocal[1] )  
     {  
         /* the case of left hand boundary domain and periodic boundary condition on its left hand boundary */  
         if( domLeft && periodic[0] ){  
             if( !periodic[2] ) {  
   
                 /* x3=0 */  
                 for( i1=0; i1<NE1; i1++ )  
                 {  
                     k = i1+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i1;  
                     node1 = numNodesLocal*N1*N2 + i1;  
   
                     out->FaceElements->Tag[k]=200;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
                     }  
                 }    
                 faceNECount  += NE1;  
                 totalNECount += NE1;  
   
                 /* x3=1 */  
                 for( i1=0; i1<NE1; i1++ )  
                 {  
                     k = i1+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i1 + (NE2-1)*N1;  
                     node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1;  
   
                     out->FaceElements->Tag[k]=200;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
                       
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1;  
                     }  
                 }    
                 faceNECount  += NE1;  
                 totalNECount += NE1;  
             }  
   
             if( !periodic[1] ) {  
                 /* x2=0 */  
                 for( i2=0; i2<NE2; i2++ )  
                 {  
                     k = i2+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i2*N1;  
                     node1 = numNodesLocal*N1*N2 + i2*N1;  
   
                     out->FaceElements->Tag[k]=20;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1+1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+N1;  
                     }  
                 }    
                 faceNECount  += NE2;  
                 totalNECount += NE2;  
   
                 /* x2=1 */  
                 for( i2=0; i2<NE2; i2++ )  
                 {  
                     k = i2+faceNECount;  
                     node0 = numDOFLocal*N1*N2 + i2*N1 + N1-2;  
                     node1 = numNodesLocal*N1*N2 + i2*N1 + N1-2;  
   
                     out->FaceElements->Tag[k]=20;  
                     out->FaceElements->Id[k]=idCount++;  
                     out->FaceElements->Color[k]=0;  
   
                     if( useElementsOnFace )  
                     {  
                         out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node1;  
                         out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1;  
                         out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+N1;  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                     } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node1+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
                         out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1+1;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1+1;  
                     }  
                 }    
                 faceNECount  += NE2;  
                 totalNECount += NE2;  
             }  
               
         }  
         if( !periodic[2] ) {  
             /* x3=0 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numDOFLocal*(i1+1) - 1;  
                 node1 = (numNodesLocal+periodicLocal[0]+periodicLocal[1])*N1*N2 + i1 + domInternal*N1*N2;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1+N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
   
             /* x3=1 */  
             for( i1=0; i1<NE1; i1++ )  
             {  
                 k = i1+faceNECount;  
                 node0 = numDOFLocal*N1*(NE2-1) + (i1+1)*numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i1 + (NE2-1)*N1 + (domInternal+periodicLocal[1]+periodicLocal[0])*N1*N2;  
   
                 out->FaceElements->Tag[k]=200;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
                   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal*N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numDOFLocal*N1+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1+N1;  
                 }  
             }    
             faceNECount  += NE1;  
             totalNECount += NE1;  
         }  
         if( !periodic[1] ) {  
             /* x2=0 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = N1*numDOFLocal*i2 + numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
   
                 if( useElementsOnFace )  
                 {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*numDOFLocal;  
                 }  
             }    
             faceNECount  += NE2;  
             totalNECount += NE2;  
   
             /* x2=1 */  
             for( i2=0; i2<NE2; i2++ )  
             {  
                 k = i2+faceNECount;  
                 node0 = numDOFLocal*N1*i2 + NE1*numDOFLocal - 1;  
                 node1 = numNodesLocal*N1*N2 + i2*N1 + (NE1-1) + (domInternal+periodicLocal[0]+periodicLocal[1])*N1*N2;  
   
                 out->FaceElements->Tag[k]=20;  
                 out->FaceElements->Id[k]=idCount++;  
                 out->FaceElements->Color[k]=0;  
                   
                 if( useElementsOnFace ){  
                     out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node1;  
                     out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node1+N1;  
                     out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N1*numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numDOFLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node1+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node1+N1+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N1*numDOFLocal+numDOFLocal;  
                 }  
             }  
             faceNECount  += NE2;  
             totalNECount += NE2;  
         }  
     }  
1020    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
1021    out->FaceElements->maxColor=0;//23;    out->FaceElements->maxColor=23;
1022        out->FaceElements->numElements=faceNECount;
1023        
1024        Finley_ElementFile_setDomainFlags( out->FaceElements );
1025    
1026      out->FaceElements->elementDistribution->numBoundary = faceNECount - out->FaceElements->elementDistribution->numInternal;    /* setup distribution info for other elements */
1027    out->FaceElements->elementDistribution->numLocal = faceNECount;      Finley_ElementFile_setDomainFlags( out->ContactElements );
1028        Finley_ElementFile_setDomainFlags( out->Points );
1029    
1030        /* reorder the degrees of freedom */
1031        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1032    
   /* setup distribution info for other elements */  
   out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;  
   out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;  
       
1033    /*   condense the nodes: */    /*   condense the nodes: */
1034    Finley_Mesh_resolveNodeIds( out );    Finley_Mesh_resolveNodeIds(out);
1035      if( !Finley_MPI_noError(mpi_info) )
1036    /* setup the CommBuffer */    {
   Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
   if ( !Finley_MPI_noError( mpi_info )) {  
     if( Finley_noError() )  
       Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );  
1037      Paso_MPIInfo_dealloc( mpi_info );      Paso_MPIInfo_dealloc( mpi_info );
1038      Finley_Mesh_dealloc(out);      Finley_Mesh_dealloc(out);
1039      return NULL;      return NULL;
1040    }    }
   
   Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
1041    
1042    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
1043    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out);
1044      if( !Finley_MPI_noError(mpi_info) )
1045      {
1046        Paso_MPIInfo_dealloc( mpi_info );
1047        Finley_Mesh_dealloc(out);
1048        return NULL;
1049      }
1050    
1051  //  print_mesh_statistics( out );    /* free up memory */
1052      Paso_MPIInfo_dealloc( mpi_info );
1053    
1054    #ifdef Finley_TRACE    #ifdef Finley_TRACE
1055    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1056    #endif    #endif
       
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
1057    
1058    return out;    return out;
1059  }  }
1060  #endif  #endif
1061    

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

  ViewVC Help
Powered by ViewVC 1.1.26