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

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

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

revision 1222 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 26  Line 26 
26    
27  #include "RectangularMesh.h"  #include "RectangularMesh.h"
28    
 #ifdef PASO_MPI  
 /* get the number of nodes/elements for domain with rank=rank, of size processors  
    where n is the total number of nodes/elements in the global domain */  
 static index_t domain_MODdim( index_t rank, index_t size, index_t n )  
 {  
   rank = size-rank-1;  
   
   if( rank < n%size )  
     return (index_t)floor(n/size)+1;  
   return (index_t)floor(n/size);  
 }  
   
   
 /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */  
 /* A bit messy, but it only has to be done once... */  
 static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal, dim_t *DOFBoundary )  
 {  
   index_t i0;  
   dim_t numNodesGlobal = numElementsGlobal+1;  
   
   (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;  
   if( rank<size-1 ) // add on node for right hand boundary  
     (*numNodesLocal) += 1;  
   
   numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;  
   periodicLocal[0] = periodicLocal[1] = FALSE;  
   nodesExternal[0] = nodesExternal[1] = 1;  
   if( periodic )  
   {  
     if( size==1 )  
     {  
       numElementsLocal[0] = numElementsGlobal;  
       nodesExternal[0] = nodesExternal[1] = 0;  
       periodicLocal[0] = periodicLocal[1] = TRUE;  
     }  
     else  
     {  
       if( rank==0 )  
       {  
         periodicLocal[0] = TRUE;  
         numNodesLocal[0]++;  
       }  
       if( rank==(size-1) )  
       {        
         periodicLocal[1] = TRUE;  
         numNodesLocal[0]--;  
         numElementsLocal[0]--;  
       }  
     }  
   }  
   else if( !periodic )  
   {  
     if( rank==0 ){  
       nodesExternal[0]--;  
       numElementsLocal[0]--;  
     }  
     if( rank==(size-1) )  
     {  
       nodesExternal[1]--;  
       numElementsLocal[0]--;  
     }  
   }  
   nodesExternal[0]*=2;  
   numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];  
     
   numElementsInternal[0] = numElementsLocal[0];  
   if( (rank==0) && (rank==size-1) );  
       
   else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )  
       numElementsInternal[0] -= 1;  
   else  
     numElementsInternal[0] -= 2;  
   
   firstNode[0] = 0;  
   for( i0=0; i0<rank; i0++ )  
     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );  
   firstNode[0] *= 2;  
     
   numDOFLocal[0] = numNodesLocal[0];  
   if( periodicLocal[0] )  
     numDOFLocal[0]--;  
     
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
   
   if( size>1 )  
   {  
     DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;  
     DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;  
   }  
   else  
   {  
     DOFBoundary[0] = DOFBoundary[1] = 0;  
   }  
 }  
   
   
 #endif  
 /**************************************************************/  
29    
30  #ifdef PASO_MPI  Finley_Mesh* Finley_RectangularMesh_Rec8(dim_t* numElements,
31  Finley_Mesh* Finley_RectangularMesh_Rec8_singleCPU(int* numElements,double* Length,int* periodic,int order, index_t reduced_order, int useElementsOnFace, Paso_MPIInfo *mpi_info)                                            double* Length,
32  #else                                            bool_t* periodic,
33  Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order, index_t reduced_order, int useElementsOnFace)                                            index_t order,
34  #endif                                            index_t reduced_order,
35                                              bool_t useElementsOnFace,
36                                              bool_t useFullElementOrder)
37  {  {
38    dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;    #define N_PER_E 2
39    index_t k,node0;    #define DIM 2
40    Finley_Mesh* out=NULL;    dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1;
41      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN, local_NE0, local_NE1, local_N0, local_N1;
42      index_t e_offset1, e_offset0, offset0, offset1, global_i0, global_i1;
43      index_t node0, myRank;
44      Finley_Mesh* out;
45      Paso_MPIInfo *mpi_info = NULL;
46    char name[50];    char name[50];
47    double time0=Finley_timer();    double time0=Finley_timer();
   NE0=MAX(1,numElements[0]);  
   NE1=MAX(1,numElements[1]);  
   N0=2*NE0+1;  
   N1=2*NE1+1;  
48    
49  #ifndef PASO_MPI    /* get MPI information */
50    if (N0<=N1) {    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
51       M0=1;    if (! Finley_noError()) {
52       M1=N0;          return NULL;
   } else {  
      M0=N1;  
      M1=1;  
   }  
 #else  
      M0=1;  
      M1=N0;  
 #endif  
   
   NFaceElements=0;  
   if (!periodic[0]) {  
       NDOF0=N0;  
       NFaceElements+=2*NE1;  
   } else {  
       NDOF0=N0-1;  
   }  
   if (!periodic[1]) {  
       NDOF1=N1;  
       NFaceElements+=2*NE0;  
   } else {  
       NDOF1=N1-1;  
53    }    }
54      myRank=mpi_info->rank;
55    
56    /*  allocate mesh: */    /*  allocate mesh: */  
     
57    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
58  #ifdef PASO_MPI    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
59    out=Finley_Mesh_alloc(name,2,order, reduced_order,mpi_info);    if (! Finley_noError()) {
60  #else        Paso_MPIInfo_free( mpi_info );
   out=Finley_Mesh_alloc(name,2,order, reduced_order);  
 #endif  
   if (! Finley_noError()) return NULL;  
   
 #ifdef PASO_MPI  
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order,mpi_info);  
 #else  
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order);  
 #endif  
   
   if (! Finley_noError()) {  
        Finley_Mesh_dealloc(out);  
        return NULL;  
   }  
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0*N1);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
 #ifdef PASO_MPI  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1, 0, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1, NE0*NE1);  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
61        return NULL;        return NULL;
62    }    }
63    
64    /*  set nodes: */    if (useFullElementOrder) {
65    #pragma omp parallel for private(i0,i1,k)       Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec9,
66    for (i1=0;i1<N1;i1++) {                                              out->order,
67      for (i0=0;i0<N0;i0++) {                                              out->reduced_order,
68        k=M0*i0+M1*i1;                                              mpi_info));
69        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];       if (useElementsOnFace) {
70        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];           Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
71        out->Nodes->Id[k]=i0+N0*i1;       } else {
72        out->Nodes->Tag[k]=0;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
73        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);                                                      out->order,
74  #ifdef PASO_MPI                                                      out->reduced_order,
75        out->Nodes->Dom[k]=NODE_INTERNAL;                                                      mpi_info));
76  #endif           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
77      }                                                         out->order,
78    }                                                         out->reduced_order,
79    /* tags for the faces: */                                                         mpi_info));
   if (!periodic[1]) {  
     for (i0=0;i0<N0;i0++) {  
       out->Nodes->Tag[M0*i0+M1*0]+=10;  
       out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;  
     }  
   }  
   if (!periodic[0]) {  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[M0*0+M1*i1]+=1;  
       out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;  
     }  
   }  
   
   /*   set the elements: */  
     
   #pragma omp parallel for private(i0,i1,k,node0)  
   for (i1=0;i1<NE1;i1++) {  
     for (i0=0;i0<NE0;i0++) {  
       k=i0+NE0*i1;  
       node0=2*i0+2*i1*N0;  
   
       out->Elements->Id[k]=k;  
       out->Elements->Tag[k]=0;  
       out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);  
 #ifdef PASO_MPI  
       out->Elements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
       out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;  
       out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;  
       out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;  
       out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;  
       out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;  
       out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;  
       out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;  
   
     }  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     
   /*   face elements: */  
   if (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=3;  
   }  
     
   totalNECount=NE0*NE1;  
   faceNECount=0;  
     
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x2=0): */  
     
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
        k=i0+faceNECount;  
        node0=2*i0;  
   
        out->FaceElements->Id[k]=i0+totalNECount;  
        out->FaceElements->Tag[k]=10;  
        out->FaceElements->Color[k]=i0%2;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
        }  
80       }       }
81       totalNECount+=NE0;    } else  {
82       faceNECount+=NE0;       Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info));
83           if (useElementsOnFace) {
84       /* **  elements on boundary 020 (x2=1): */           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8Face,
85                                                                        out->order,
86       #pragma omp parallel for private(i0,k,node0)                                                                    out->reduced_order,
87       for (i0=0;i0<NE0;i0++) {                                                                    mpi_info));
88         k=i0+faceNECount;           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8Face_Contact,
89         node0=2*i0+2*(NE1-1)*N0;                                                                      out->order,
90                                                                        out->reduced_order,
91         out->FaceElements->Id[k]=i0+totalNECount;                                                                      mpi_info));
92         out->FaceElements->Tag[k]=20;       } else {
93         out->FaceElements->Color[k]=i0%2+2;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
94  #ifdef PASO_MPI                                                                    out->order,
95           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                                                                    out->reduced_order,
96  #endif                                                                    mpi_info));
97         if (useElementsOnFace) {           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
98            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;                                                                       out->order,
99            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;                                                                       out->reduced_order,
100            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;                                                                       mpi_info));
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;  
        }  
101       }       }
      totalNECount+=NE0;  
      faceNECount+=NE0;  
102    }    }
103    if (!periodic[0]) {    Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
104       /* **  elements on boundary 001 (x1=0): */                                                   out->order,
105                                                       out->reduced_order,
106       #pragma omp parallel for private(i1,k,node0)                                                   mpi_info));
      for (i1=0;i1<NE1;i1++) {  
        k=i1+faceNECount;  
        node0=2*i1*N0;  
   
        out->FaceElements->Id[k]=i1+totalNECount;  
        out->FaceElements->Tag[k]=1;  
        out->FaceElements->Color[k]=(i1%2)+4;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;  
        }  
     
      }  
      totalNECount+=NE1;  
      faceNECount+=NE1;  
     
      /* **  elements on boundary 002 (x1=1): */  
       
      #pragma omp parallel for private(i1,k,node0)  
      for (i1=0;i1<NE1;i1++) {  
        k=i1+faceNECount;  
        node0=2*(NE0-1)+2*i1*N0;  
     
        out->FaceElements->Id[k]=i1+totalNECount;  
        out->FaceElements->Tag[k]=2;  
        out->FaceElements->Color[k]=(i1%2)+6;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
     
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;  
        }  
      }  
      totalNECount+=NE1;  
      faceNECount+=NE1;  
   
   }  
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=7;  
   
 #ifdef PASO_MPI  
     Finley_ElementFile_setDomainFlags( out->Elements );  
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     /* reorder the degrees of freedom */  
     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );  
 #endif  
   
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
       
107    if (! Finley_noError()) {    if (! Finley_noError()) {
108        Finley_Mesh_dealloc(out);        Paso_MPIInfo_free( mpi_info );
109          Finley_Mesh_free(out);
110        return NULL;        return NULL;
111    }    }
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
   
   
   #ifdef Finley_TRACE  
   printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);  
   #endif  
112    
113    if (Finley_noError()) {    /* set up the global dimensions of the mesh */
        if ( ! Finley_Mesh_isPrepared(out)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
       Finley_Mesh_dealloc(out);  
   }  
   return out;  
 }  
   
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order, index_t reduced_order ,int useElementsOnFace) {  
   dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];  
     dim_t N0t, NDOF0t;  
   index_t *numForward=NULL, *numBackward=NULL;  
   index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, targetDomain=-1,faceNEBoundary;  
   Finley_Mesh* out=NULL;  
   char name[50];  
     index_t *indexForward=NULL, *indexBackward=NULL;  
   double time0=Finley_timer();  
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;  
   Paso_MPIInfo *mpi_info = NULL;  
     index_t faceStart, numFaceLeft;  
   
     /* these are used in constructing the face elements, and give the permutation  
        of the node order of a normal element for each face */  
     index_t *facePerm;    
     index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};  
     index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};  
     index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};  
     index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};  
     index_t face0nc[]={0, 1, 4};  
     index_t face1nc[]={2, 3, 6};  
     index_t face2nc[]={3, 0, 7};  
     index_t face3nc[]={1, 2, 5};  
114    
115    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
116    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
117    N0=2*NE0+1;    N0=N_PER_E*NE0+1;
118    N1=2*NE1+1;    N1=N_PER_E*NE1+1;
119    
120    /* get MPI information */    /* work out the largest dimension */
121    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    if (N1=MAX(N0,N1)) {
122    if (! Finley_noError()) {       Nstride0=1;
123          Finley_Mesh_dealloc(out);       Nstride1=N0;
124          return NULL;       local_NE0=NE0;
125    }       e_offset0=0;
126         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
127      } else {
128         Nstride0=N1;
129         Nstride1=1;
130         Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
131         local_NE1=NE1;
132         e_offset1=0;
133      }
134      offset0=e_offset0*N_PER_E;
135      offset1=e_offset1*N_PER_E;
136      local_N0=local_NE0*N_PER_E+1;
137      local_N1=local_NE1*N_PER_E+1;
138    
139      // use the serial method for the single CPU case    /* get the number of surface elements */
     if( mpi_info->size==1 ){  
         out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order, reduced_order,useElementsOnFace, mpi_info);  
         return out;  
     }  
       
   if( mpi_info->rank==0 )  
     domLeft = TRUE;  
   if( mpi_info->rank==mpi_info->size-1 )  
     domRight = TRUE;  
   if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )  
     domInternal = TRUE;  
   
   /* dimensions of the local subdomain */  
   domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );    
140    
141    NFaceElements=0;    NFaceElements=0;
142    if (!periodic[0]) {    if (!periodic[0]) {
143        NDOF0=N0;       NDOF0=N0;
144        NFaceElements+=(domLeft + domRight)*NE1;       if (e_offset0 == 0) NFaceElements+=local_NE1;
145         if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
146    } else {    } else {
147        NDOF0=N0-1;        NDOF0=N0-1;
148    }    }
149    if (!periodic[1]) {    if (!periodic[1]) {
150        NDOF1=N1;       NDOF1=N1;
151        NFaceElements+=2*numElementsLocal;       if (e_offset1 == 0) NFaceElements+=local_NE0;
152         if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
153    } else {    } else {
154        NDOF1=N1-1;        NDOF1=N1-1;
155    }    }
156    
     boundaryLeft = !domLeft || periodicLocal[0];  
     boundaryRight = !domRight || periodicLocal[1];  
     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;  
     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;  
     firstNodeConstruct = firstNode - 2*boundaryLeft;  
     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;  
   
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d mesh",N0,N1);  
   out=Finley_Mesh_alloc(name,2,order,reduced_order,mpi_info);  
   
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order,mpi_info);  
   
   if (! Finley_noError()) {  
        Finley_Mesh_dealloc(out);  
        return NULL;  
   }  
     
157    /*  allocate tables: */    /*  allocate tables: */
158    Finley_NodeFile_allocTable(out->Nodes,N0t*N1);  
159    Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
160      Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
161    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
       
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1, NE1*(numElementsLocal-boundaryRight*(!periodic[1])) );  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(!periodic[1]) );  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
162    
163    /*  set nodes: */    if (Finley_noError()) {
164    #pragma omp parallel for private(i0,i1,k)       /* create nodes */
165    for (i1=0;i1<N1;i1++) {    
166      for (i0=0;i0<N0t;i0++, k++) {       #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
167              k=i0+i1*N0t;       for (i1=0;i1<local_N1;i1++) {
168        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];         for (i0=0;i0<local_N0;i0++) {
169        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];             k=i0+local_N0*i1;
170        out->Nodes->Id[k]=i0 + i1*N0t;             global_i0=i0+offset0;
171        out->Nodes->Tag[k]=0;             global_i1=i1+offset1;
172        out->Nodes->Dom[k]=NODE_INTERNAL;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
173        out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
174      }             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
175          /* take care to make the DOF reflect the internal/external position of the DOF */             out->Nodes->Tag[k]=0;
176          if( boundaryLeft ){             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
177              out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;                                                     +Nstride1*(global_i1%NDOF1);
178              out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;           }
             out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;    
             if( periodicLocal[0] ){  
                 out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];  
                 out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;    
             }  
         }  
         if( boundaryRight ){  
             out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;    
             out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;    
             out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;    
         }  
   }  
   /* tags for the faces: */  
   if (!periodic[1]) {  
     for (i0=0;i0<N0t;i0++) {  
       out->Nodes->Tag[i0]+=10;  
       out->Nodes->Tag[i0+N0t*(N1-1)]+=20;  
     }  
   }  
   if (domLeft && !periodic[0]) {  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[i1*N0t]+=1;  
     }  
   }  
     if (domRight && !periodic[0]){  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[(i1+1)*N0t-1]+=2;  
         }  
     }  
     
     /* setup the receive information */  
     indexBackward = TMPMEMALLOC( NDOF1*2, index_t );  
     indexForward  = TMPMEMALLOC( NDOF1*2, index_t );  
     if( !(mpi_info->size==2 && periodicLocal[0])){  
         /* Backward (receive) DOF */  
         if( boundaryLeft ){  
             targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;  
   
             for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];  
                 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];  
                 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];  
             }  
   
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );  
         }  
           
         if( boundaryRight ){  
             targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;  
   
             for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];  
                 indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];  
                 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];  
             }  
   
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );  
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );  
         }  
     }  
     else{  
         targetDomain = 1;  
   
         for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];  
                 indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];  
                 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];  
         }  
   
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );  
   
         for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];  
                 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];  
                 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];  
         }  
   
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );  
     }  
     TMPMEMFREE( indexBackward );  
     TMPMEMFREE( indexForward );  
   
   /*   set the elements: */  
   #pragma omp parallel for private(i0,i1,k,node0)  
   for (i1=0;i1<NE1;i1++) {  
     for (i0=0;i0<numElementsLocal;i0++) {  
       k=i0+numElementsLocal*i1;  
             node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t :  2*i0+2*i1*N0t+periodicLocal[0] ;  
   
       out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0) * NE1 + i1;  
       out->Elements->Tag[k]=0;  
       out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);  
             out->Elements->Dom[k] = ELEMENT_INTERNAL;  
   
       out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;  
       out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;  
       out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;  
       out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;  
       out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;  
       out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;  
       out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;  
     }  
         if( boundaryLeft )  
             out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;  
         if( boundaryRight )  
             out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     
     /* set the distribution information for the elements */  
     Finley_ElementFile_setDomainFlags( out->Elements );  
   
   /*   face elements: */  
   if (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=3;  
   }  
     
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
   totalNECount=numElementsLocal*NE1;  
   faceNECount=0;  
     faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);  
     
     numFaceLeft = domLeft*(!periodic[0])*NE1;  
     faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];  
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x1=0): */  
     
      #pragma omp parallel for private(i0,k,kk,facePerm)  
      for (i0=0;i0<numElementsLocal; i0++) {  
        k=i0+faceNECount;  
              kk=i0;  
              facePerm = useElementsOnFace ? face0 : face0nc;  
   
        out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;  
        out->FaceElements->Tag[k]=10;  
        out->FaceElements->Color[k]=i0%2;  
              out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
                   
              for( j=0; j<NUMNODES; j++ )  
                   out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
179       }       }
180           if( boundaryLeft ){       /*   set the elements: */
181               out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;       NN=out->Elements->numNodes;
182               if( periodicLocal[0] )       #pragma omp parallel for private(i0,i1,k,node0)
183                   out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;       for (i1=0;i1<local_NE1;i1++) {
184           }           for (i0=0;i0<local_NE0;i0++) {
185           if( boundaryRight ){            
186               out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;             k=i0+local_NE0*i1;        
187               out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
188           }  
189       totalNECount+=numElementsLocal;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
190       faceNECount+=numElementsLocal;             out->Elements->Tag[k]=0;
191                 out->Elements->Owner[k]=myRank;
192       /* **  elements on boundary 020 (x1=1): */  
193                 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
194       #pragma omp parallel for private(i0,k,kk,facePerm)             out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
195       for (i0=0;i0<numElementsLocal;i0++) {             out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
196         k=i0+faceNECount;             out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
197               kk=i0 + numElementsLocal*(NE1-1);             out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
198               facePerm = useElementsOnFace ? face1 : face1nc;             out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
199               out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
200         out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;             out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
201         out->FaceElements->Tag[k]=20;             if (useFullElementOrder) {
202         out->FaceElements->Color[k]=i0%2+2;                out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
203               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             }
204             }
              for( j=0; j<NUMNODES; j++ )  
                   out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
205       }       }
206           if( boundaryLeft ){       /* face elements */
207               out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;       NN=out->FaceElements->numNodes;
208               if( periodicLocal[0] )       totalNECount=NE0*NE1;
209                   out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;       faceNECount=0;
210           }       if (!periodic[0]) {
211           if( boundaryRight ){          /* **  elements on boundary 001 (x1=0): */
212               out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;      
213               out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;          if (e_offset0 == 0) {
214           }             #pragma omp parallel for private(i1,k,node0)
215       totalNECount+=numElementsLocal;             for (i1=0;i1<local_NE1;i1++) {
216       faceNECount+=numElementsLocal;        
217    }                 k=i1+faceNECount;
218    if (domLeft && !periodicLocal[0]) {                 node0=Nstride1*N_PER_E*(i1+e_offset1);
219       /* **  elements on boundary 001 (x0=0): */  
220                     out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
221       #pragma omp parallel for private(i1,k,kk,facePerm)                 out->FaceElements->Tag[k]=1;
222       for (i1=0;i1<NE1;i1++) {                 out->FaceElements->Owner[k]=myRank;
223         k=i1+faceNECount;                 if (useElementsOnFace) {
224               kk=i1*numElementsLocal;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
225               facePerm = useElementsOnFace ? face2 : face2nc;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
226                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
227         out->FaceElements->Id[k]=faceStart + i1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
228         out->FaceElements->Tag[k]=1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
229         out->FaceElements->Color[k]=(i1%2)+4;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
230               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
231                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
232               for( j=0; j<NUMNODES; j++ )                 } else {
233                    out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
234                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
235                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
236                   }
237               }
238               faceNECount+=local_NE1;
239            }
240            totalNECount+=NE1;
241            /* **  elements on boundary 002 (x1=1): */
242            if (local_NE0+e_offset0 == NE0) {
243               #pragma omp parallel for private(i1,k,node0)
244               for (i1=0;i1<local_NE1;i1++) {
245                   k=i1+faceNECount;
246                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
247    
248                   out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
249                   out->FaceElements->Tag[k]=2;
250                   out->FaceElements->Owner[k]=myRank;
251    
252                   if (useElementsOnFace) {
253                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
254                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
255                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
256                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
257                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
258                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
259                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
260                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
261                   } else {
262                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
263                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
264                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
265                   }
266               }
267               faceNECount+=local_NE1;
268             }
269             totalNECount+=NE1;
270       }       }
271       totalNECount+=NE1;       if (!periodic[1]) {
272       faceNECount+=NE1;          /* **  elements on boundary 010 (x2=0): */
273              if (e_offset1 == 0) {
274       /* **  elements on boundary 002 (x0=1): */             #pragma omp parallel for private(i0,k,node0)
275    }             for (i0=0;i0<local_NE0;i0++) {
276    if (domRight && !periodicLocal[1]) {                 k=i0+faceNECount;
277       #pragma omp parallel for private(i1,k,kk,facePerm)                 node0=Nstride0*N_PER_E*(i0+e_offset0);
278       for (i1=0;i1<NE1;i1++) {          
279         k=i1+faceNECount;                 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
280               kk=(i1+1)*numElementsLocal-1;                 out->FaceElements->Tag[k]=10;
281               facePerm = useElementsOnFace ? face3 : face3nc;                 out->FaceElements->Owner[k]=myRank;
282          
283                   if (useElementsOnFace) {
284                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
285                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
286                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
287                       out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
288                       out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
289                       out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
290                       out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
291                       out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
292                   } else {
293                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
294                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
295                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
296                   }
297               }
298               faceNECount+=local_NE0;
299            }
300            totalNECount+=NE0;
301            /* **  elements on boundary 020 (x2=1): */
302            if (local_NE1+e_offset1 == NE1) {
303               #pragma omp parallel for private(i0,k,node0)
304               for (i0=0;i0<local_NE0;i0++) {
305                   k=i0+faceNECount;
306                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
307        
308         out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;                 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
309         out->FaceElements->Tag[k]=2;                 out->FaceElements->Tag[k]=20;
310         out->FaceElements->Color[k]=(i1%2)+6;                 out->FaceElements->Owner[k]=myRank;
311               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                 if (useElementsOnFace) {
312                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
313                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
314                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
315                        out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
316                        out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
317                        out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
318                        out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
319                        out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
320                   } else {
321                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
322                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
323                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
324                   }
325               }
326               faceNECount+=local_NE0;
327            }
328            totalNECount+=NE0;
329         }
330         /* add tag names */
331         Finley_Mesh_addTagMap(out,"top", 20);
332         Finley_Mesh_addTagMap(out,"bottom", 10);
333         Finley_Mesh_addTagMap(out,"left", 1);
334         Finley_Mesh_addTagMap(out,"right", 2);
335        
336               for( j=0; j<NUMNODES; j++ )       /* prepare mesh for further calculatuions:*/
337                    out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];       if (Finley_noError()) {
338             Finley_Mesh_resolveNodeIds(out);
339         }
340         if (Finley_noError()) {
341             Finley_Mesh_prepare(out);
342       }       }
      totalNECount+=NE1;  
      faceNECount+=NE1;  
343    }    }
344    out->FaceElements->minColor=0;    if (!Finley_noError()) {
345    out->FaceElements->maxColor=7;        Finley_Mesh_free(out);
     out->FaceElements->numElements = faceNECount;  
   
   /* setup distribution info for other elements */  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     Finley_Mesh_prepareElementDistribution( out );  
   
     /* reorder the degrees of freedom */  
     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );  
       
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
   if ( !Finley_MPI_noError( mpi_info )) {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
346    }    }
347      /* free up memory */
348    /* prepare mesh for further calculatuions:*/    Paso_MPIInfo_free( mpi_info );  
   Finley_Mesh_prepare(out) ;  
   
349    #ifdef Finley_TRACE    #ifdef Finley_TRACE
350    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
351    #endif    #endif
352    
   if ( !Finley_MPI_noError( mpi_info )) {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
   if (Finley_noError()) {  
        if (! Finley_Mesh_isPrepared(out)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   }  
353    return out;    return out;
354  }  }
 #endif  

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

  ViewVC Help
Powered by ViewVC 1.1.26