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

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

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

revision 1222 by ksteube, Tue May 15 03:23:17 2007 UTC revision 1223 by gross, Fri Aug 3 02:40:39 2007 UTC
# Line 27  Line 27 
27    
28  #include "RectangularMesh.h"  #include "RectangularMesh.h"
29    
 /**************************************************************/  
30    
31  #ifdef PASO_MPI  Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,
32  /* get the number of nodes/elements for domain with rank=rank, of size processors                                            double* Length,
33     where n is the total number of nodes/elements in the global domain */                                            bool_t* periodic,
34  static index_t domain_MODdim( index_t rank, index_t size, index_t n )                                            index_t order,
35  {                                            index_t reduced_order,
36    rank = size-rank-1;                                            bool_t useElementsOnFace,
37                                              bool_t useFullElementOrder)
   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 */  
 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 *numDOFInternal )  
 {  
   index_t i0;  
   dim_t numNodesGlobal = numElementsGlobal+1;  
   
   (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );  
   
   numElementsLocal[0] = numNodesLocal[0]+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]--;  
     }  
   }  
   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 );  
   
   numDOFLocal[0] = numNodesLocal[0];  
   if( periodicLocal[0] )  
     numDOFLocal[0]--;  
           
     numDOFInternal[0] = numDOFLocal[0];  
     if( size>1 )  
     {    
         if( rank>0 || periodic ){  
             numDOFInternal[0] -= 1;  
         }    
         if( rank<(size-1) || periodic ){  
             numDOFInternal[0] -= 1;  
         }    
     }    
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
 }  
 #endif  
   
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Rec4_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)  
 #else  
 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace)  
 #endif  
38  {  {
39    dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;    #define N_PER_E 1
40    index_t k,node0;    #define DIM 2
41      dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1, local_NE0, local_NE1, local_N0, local_N1, global_i0, global_i1;
42      index_t offset0, offset1, e_offset0, e_offset1;
43      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN;
44      index_t node0, myRank;
45    Finley_Mesh* out;    Finley_Mesh* out;
46      Paso_MPIInfo *mpi_info = NULL;
47    char name[50];    char name[50];
48    double time0=Finley_timer();    double time0=Finley_timer();
   NE0=MAX(1,numElements[0]);  
   NE1=MAX(1,numElements[1]);  
   N0=NE0+1;  
   N1=NE1+1;  
49    
50    if (N0<=N1) {    /* get MPI information */
51       M0=1;    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
52       M1=N0;    if (! Finley_noError()) {
53    } else {          return NULL;
      M0=N1;  
      M1=1;  
   }  
   
   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;  
54    }    }
55      myRank=mpi_info->rank;
56    
57    /*  allocate mesh: */    /*  allocate mesh: */  
     
58    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
59  #ifdef PASO_MPI    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
60    out=Finley_Mesh_alloc(name,2,order,reduced_order,mpi_info);    if (! Finley_noError()) {
61          Paso_MPIInfo_free( mpi_info );
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Rec4,out->order, out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1, out->reduced_order,out->order,mpi_info);  
   
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
62        return NULL;        return NULL;
63    }    }
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0*N1);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   
   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 );  
 #else  
   out=Finley_Mesh_alloc(name,2,order,reduced_order);  
64    
65    if (! Finley_noError()) return NULL;    Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec4,out->order,out->reduced_order,mpi_info));
   
   out->Elements=Finley_ElementFile_alloc(Rec4,out->order, out->reduced_order);  
66    if (useElementsOnFace) {    if (useElementsOnFace) {
67       out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, out->reduced_order);           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec4Face,
68       out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, out->reduced_order);                                                                    out->order,
69    } else {                                                                    out->reduced_order,
70       out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, out->reduced_order);                                                                    mpi_info));
71       out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, out->reduced_order);           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec4Face_Contact,
72    }                                                                      out->order,
73    out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order);                                                                      out->reduced_order,
74                                                                        mpi_info));
75    if (! Finley_noError()) {    } else {
76        Finley_Mesh_dealloc(out);           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line2,
77        return NULL;                                                                    out->order,
78    }                                                                    out->reduced_order,
79                                                                        mpi_info));
80    /*  allocate tables: */           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line2_Contact,
81    Finley_NodeFile_allocTable(out->Nodes,N0*N1);                                                                       out->order,
82    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);                                                                       out->reduced_order,
83    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);                                                                       mpi_info));
84  #endif    }
85      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
86                                                     out->order,
87                                                     out->reduced_order,
88                                                     mpi_info));
89    if (! Finley_noError()) {    if (! Finley_noError()) {
90        Finley_Mesh_dealloc(out);        Paso_MPIInfo_free( mpi_info );
91          Finley_Mesh_free(out);
92        return NULL;        return NULL;
93    }    }
   /*  set nodes: */  
   #pragma omp parallel for private(i0,i1,k)  
   for (i1=0;i1<N1;i1++) {  
     for (i0=0;i0<N0;i0++) {  
       k=M0*i0+M1*i1;  
       out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];  
       out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
       out->Nodes->Id[k]=i0+N0*i1;  
       out->Nodes->Tag[k]=0;  
       out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);  
 #ifdef PASO_MPI  
       out->Nodes->Dom[k]=NODE_INTERNAL;  
 #endif  
     }  
   }  
   /* tags for the faces: */  
   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=i0+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,4)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;  
       out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0+1;  
       out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0;  
   
     }  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     
   /*   face elements: */  
     
   if (useElementsOnFace) {  
      NUMNODES=4;  
   } else {  
      NUMNODES=2;  
   }  
   totalNECount=NE0*NE1;  
   faceNECount=0;  
     
   /* elements on boundary 010 (x2=0): */  
   if (!periodic[1]) {  
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
           k=i0+faceNECount;  
           node0=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+1;  
               out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;  
               out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;  
           } else {  
               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
               out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
           }  
      }  
      totalNECount+=NE0;  
      faceNECount+=NE0;  
       
      /* **  elements on boundary 020 (x2=1): */  
     
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
          k=i0+faceNECount;  
          node0=i0+(NE1-1)*N0;  
     
          out->FaceElements->Id[k]=i0+totalNECount;  
          out->FaceElements->Tag[k]=20;  
          out->FaceElements->Color[k]=i0%2+2;  
 #ifdef PASO_MPI  
               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
     
          if (useElementsOnFace) {  
              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;  
              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;  
              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
          } else {  
              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;  
          }  
      }  
      totalNECount+=NE0;  
      faceNECount+=NE0;  
   }  
   if (!periodic[0]) {  
      /* **  elements on boundary 001 (x1=0): */  
      #pragma omp parallel for private(i1,k,node0)  
      for (i1=0;i1<NE1;i1++) {  
           k=i1+faceNECount;  
           node0=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+N0;  
               out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
               out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
               out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;  
           } else {  
               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;  
               out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
     
           }  
      }  
      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=(NE0-1)+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+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;  
              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;  
              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
           } else {  
              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;  
           }  
      }  
      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);  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
94    
95    #ifdef Finley_TRACE    /* set up the global dimensions of the mesh */
   printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);  
   #endif  
   if (Finley_noError()) {  
        if ( ! Finley_Mesh_isPrepared(out) ) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
       Finley_Mesh_dealloc(out);  
   }  
   return out;  
 }  
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order, index_t reduced_order, bool_t useElementsOnFace)  
 {  
   dim_t N0,N1,NE0,NE1,i0,i1,N0t, NDOF0t, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1,numDOFInternal;  
   index_t NUMNODES,k,firstNode=0, DOFcount=0, node0, node1;  
   index_t targetDomain=-1, firstNodeConstruct;  
   index_t *forwardDOF=NULL, *backwardDOF=NULL;  
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;  
   Finley_Mesh* out;  
   char name[50];  
   double time0=Finley_timer();  
     index_t faceStart, numFaceLeft, i;  
96    
97    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
98    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
99    N0=NE0+1;    N0=N_PER_E*NE0+1;
100    N1=NE1+1;    N1=N_PER_E*NE1+1;
   Paso_MPIInfo *mpi_info = NULL;  
101    
102    /* get MPI information */    /* work out the largest dimension */
103    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );    if (N1=MAX(N0,N1)) {
104    if (! Finley_noError())       Nstride0=1;
105          return NULL;       Nstride1=N0;
106         local_NE0=NE0;
107         e_offset0=0;
108         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
109      } else {
110         Nstride0=N1;
111         Nstride1=1;
112         Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
113         local_NE1=NE1;
114         e_offset1=0;
115      }
116      offset0=e_offset0*N_PER_E;
117      offset1=e_offset1*N_PER_E;
118      local_N0=local_NE0*N_PER_E+1;
119      local_N1=local_NE1*N_PER_E+1;
120    
121      /* use the serial code to generate the mesh in the 1-CPU case */    /* get the number of surface elements */
     if( mpi_info->size==1 ){  
         out = Finley_RectangularMesh_Rec4_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, &numDOFInternal );    
   
   if (N0<=N1) {  
      M0=1;  
      M1=N0;  
   } else {  
      M0=N1;  
      M1=1;  
   }  
122    
123    NFaceElements=0;    NFaceElements=0;
124    if (!periodic[0])    if (!periodic[0]) {
125    {       NDOF0=N0;
126        if( domLeft )       if (e_offset0 == 0) NFaceElements+=local_NE1;
127          NFaceElements+=NE1;       if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
128        if( domRight )    } else {
129          NFaceElements+=NE1;        NDOF0=N0-1;
130    }    }
131    if (!periodic[1]) {    if (!periodic[1]) {
132        NDOF1=N1;       NDOF1=N1;
133        NFaceElements+=2*numElementsLocal;       if (e_offset1 == 0) NFaceElements+=local_NE0;
134         if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
135    } else {    } else {
136        NDOF1=N1-1;        NDOF1=N1-1;
137    }    }
     
     boundaryLeft = !domLeft || periodicLocal[0];  
     boundaryRight = !domRight || periodicLocal[1];  
     N0t = numNodesLocal + boundaryRight + boundaryLeft;  
     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;  
     firstNodeConstruct = firstNode - boundaryLeft;  
     firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;  
138    
   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( Rec4, out->order, out->reduced_order, mpi_info );  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, out->reduced_order, mpi_info );  
      out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, out->reduced_order, mpi_info );  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, out->reduced_order, mpi_info );  
      out->ContactElements=Finley_ElementFile_alloc(Line2_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;  
   }  
     
139    /*  allocate tables: */    /*  allocate tables: */
   Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );  
   Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
   
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, 2*NDOF1, 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]) );  
140    
141    if (! Finley_noError()) {    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
142        Finley_Mesh_dealloc(out);    Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
143        return NULL;    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
   }  
144    
145    /*  set nodes: */    if (Finley_noError()) {
146    #pragma omp parallel for private(i0,i1,k)       /* create nodes */
147    for (i1=0;i1<N1;i1++) {       #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
148      for ( i0=0;i0<N0t;i0++,k++) {       for (i1=0;i1<local_N1;i1++) {
149              k=i0+i1*N0t;         for (i0=0;i0<local_N0;i0++) {
150        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];             k=i0+local_N0*i1;
151        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];             global_i0=i0+offset0;
152        out->Nodes->Id[k]=k;             global_i1=i1+offset1;
153        out->Nodes->Tag[k]=0;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
154        out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
155              out->Nodes->Dom[k]=NODE_INTERNAL;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
156      }             out->Nodes->Tag[k]=0;
157    }             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
158      if( boundaryLeft ){                                                 +Nstride1*(global_i1%NDOF1);
159          for( i1=0; i1<N1; i1++ ){         }
160              out->Nodes->Dom[N0t*i1] = NODE_EXTERNAL;       }
161              out->Nodes->Dom[N0t*i1+1] = NODE_BOUNDARY;       /*   set the elements: */
162          }       NN=out->Elements->numNodes;
163      }       #pragma omp parallel for private(i0,i1,k,node0)
164      if( boundaryRight ){       for (i1=0;i1<local_NE1;i1++) {
165          for( i1=0; i1<N1; i1++ ){           for (i0=0;i0<local_NE0;i0++) {
166              out->Nodes->Dom[N0t*(i1+1)-1] = NODE_EXTERNAL;            
167              out->Nodes->Dom[N0t*(i1+1)-2] = NODE_BOUNDARY;             k=i0+local_NE0*i1;        
168          }             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
169      }  
170      if( periodicLocal[0] ){             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
171          for( i1=0; i1<N1; i1++ ){             out->Elements->Tag[k]=0;
172              out->Nodes->degreeOfFreedom[i1*N0t+2] = out->Nodes->degreeOfFreedom[i1*N0t+1];             out->Elements->Owner[k]=myRank;
173              out->Nodes->Dom[N0t*i1+2] = NODE_BOUNDARY;  
174          }             out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
175      }             out->Elements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
176                       out->Elements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
177    /* tag 'em */             out->Elements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
178    for(i0=0; i0<N0t; i0++ )           }
179    {       }
180      out->Nodes->Tag[i0] += 10;                        // bottom       /* face elements */
181      out->Nodes->Tag[i0+N0t*(N1-1)] += 20;   // top       NN=out->FaceElements->numNodes;
182    }       totalNECount=NE0*NE1;
183    if( domLeft && !periodicLocal[0] )       faceNECount=0;
184      for( i1=0; i1<N1; i1++ )       if (!periodic[0]) {
185        out->Nodes->Tag[i1*N0t] += 1;         // left          /* **  elements on boundary 001 (x1=0): */
186        
187    if( domRight && !periodicLocal[0] )          if (e_offset0 == 0) {
188      for( i1=0; i1<N1; i1++ )             #pragma omp parallel for private(i1,k,node0)
189        out->Nodes->Tag[(i1+1)*N0t-1] += 2;   // right             for (i1=0;i1<local_NE1;i1++) {
190          
191      /* form the boudary communication information */                 k=i1+faceNECount;
192      forwardDOF  = MEMALLOC(NDOF1,index_t);                 node0=Nstride1*N_PER_E*(i1+e_offset1);
193      backwardDOF = MEMALLOC(NDOF1,index_t);  
194      if( !(mpi_info->size==2 && periodicLocal[0])){                 out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
195          if( boundaryLeft  ) {                 out->FaceElements->Tag[k]=1;
196              targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;                 out->FaceElements->Owner[k]=myRank;
197              for( i1=0; i1<NDOF1; i1++ ){                 if (useElementsOnFace) {
198                  forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
199                  backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
200              }                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride0;
201              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0;
202              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );                 } else {
203          }                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1;
204          if( boundaryRight ) {                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
205              targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;                 }
206              for( i1=0; i1<NDOF1; i1++ ){             }
207                  forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];             faceNECount+=local_NE1;
208                  backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];          }
209              }          totalNECount+=NE1;
210              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );          /* **  elements on boundary 002 (x1=1): */
211              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );          if (local_NE0+e_offset0 == NE0) {
212          }             #pragma omp parallel for private(i1,k,node0)
213      } else{             for (i1=0;i1<local_NE1;i1++) {
214          /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */                 k=i1+faceNECount;
215          targetDomain = 1;                 node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
216            
217          for( i1=0; i1<NDOF1; i1++ ){                 out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
218              forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];                 out->FaceElements->Tag[k]=2;
219              backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];                 out->FaceElements->Owner[k]=myRank;
220          }  
221          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );                 if (useElementsOnFace) {
222          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
223                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
224          for( i1=0; i1<NDOF1; i1++ ){                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
225              forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
226              backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];                 } else {
227          }                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride0;
228          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1+Nstride0;
229          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );                 }
230      }             }
231      MEMFREE( forwardDOF );             faceNECount+=local_NE1;
232      MEMFREE( backwardDOF );           }
233             totalNECount+=NE1;
234    /* elements: */       }
235           if (!periodic[1]) {
236   #pragma omp parallel for private(i0,i1,k,node0)            /* **  elements on boundary 010 (x2=0): */
237   for ( i1=0; i1<NE1; i1++) {          if (e_offset1 == 0) {
238      for ( i0=0; i0<numElementsLocal; i0++, k++) {             #pragma omp parallel for private(i0,k,node0)
239              k=i1*numElementsLocal + i0;             for (i0=0;i0<local_NE0;i0++) {
240              node0 = (periodicLocal[0] && !i0) ? i1*N0t :  i1*N0t + i0 + periodicLocal[0];                 k=i0+faceNECount;
241                   node0=Nstride0*N_PER_E*(i0+e_offset0);
242        out->Elements->Id[k]=((firstNodeConstruct+i0)%NE0) * NE1 + i1;          
243        out->Elements->Tag[k]=0;                 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
244        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);                 out->FaceElements->Tag[k]=10;
245              out->Elements->Dom[k]=ELEMENT_INTERNAL;                 out->FaceElements->Owner[k]=myRank;
246          
247        out->Elements->Nodes[INDEX2(0,k,4)]=node0;                 if (useElementsOnFace) {
248        out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;                     out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
249        out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0t+1;                     out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
250        out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0t;                     out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+Nstride0;
251      }                     out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1;
252                   } else {
253                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
254                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride0;
255                   }
256               }
257               faceNECount+=local_NE0;
258            }
259            totalNECount+=NE0;
260            /* **  elements on boundary 020 (x2=1): */
261            if (local_NE1+e_offset1 == NE1) {
262               #pragma omp parallel for private(i0,k,node0)
263               for (i0=0;i0<local_NE0;i0++) {
264                   k=i0+faceNECount;
265                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
266      
267                   out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
268                   out->FaceElements->Tag[k]=20;
269                   out->FaceElements->Owner[k]=myRank;
270                   if (useElementsOnFace) {
271                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
272                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
273                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
274                        out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride0;
275                   } else {
276                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+Nstride1+Nstride0;
277                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride1;
278                   }
279               }
280               faceNECount+=local_NE0;
281            }
282            totalNECount+=NE0;
283         }
284         /* add tag names */
285         Finley_Mesh_addTagMap(out,"top", 20);
286         Finley_Mesh_addTagMap(out,"bottom", 10);
287         Finley_Mesh_addTagMap(out,"left", 1);
288         Finley_Mesh_addTagMap(out,"right", 2);
289      
290         /* prepare mesh for further calculatuions:*/
291         if (Finley_noError()) {
292             Finley_Mesh_resolveNodeIds(out);
293         }
294         if (Finley_noError()) {
295             Finley_Mesh_prepare(out);
296         }
297    }    }
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     if( boundaryLeft )  
         for( i1=0; i1<NE1; i1++ )  
             out->Elements->Dom[i1*numElementsLocal]=ELEMENT_BOUNDARY;  
     if( boundaryRight )  
         for( i1=0; i1<NE1; i1++ )  
             out->Elements->Dom[(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;  
           
     Finley_ElementFile_setDomainFlags( out->Elements );  
298    
299    /*   face elements: */    if (!Finley_noError()) {
300    if (useElementsOnFace) {        Finley_Mesh_free(out);
      NUMNODES=4;  
   } else {  
      NUMNODES=2;  
301    }    }
   totalNECount=numElementsLocal*NE1;  
   faceNECount=0;  
     
     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,node0)  
       for (i0=0;i0<numElementsLocal;i0++) {  
                 k=i0+faceNECount;  
                 node0 = (periodicLocal[0] && !i0) ? 0 :  i0 + periodicLocal[0];  
   
                 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;  
   
                 if (useElementsOnFace) {  
                         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+N0t+1;  
                         out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t;  
                 } else {  
                         out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                         out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;  
                 }  
       }  
       if( boundaryLeft ){  
           out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;  
           if( periodicLocal[0] )  
               out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;  
         }  
       if( boundaryRight ){  
           out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;  
             out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];  
         }  
   
       totalNECount += numElementsLocal;  
       faceNECount  += numElementsLocal;  
       
       /* **  elements on boundary 020 (x1=1): */  
   
       #pragma omp parallel for private(i0,k,node0)  
       for (i0=0;i0<numElementsLocal;i0++) {  
             k=i0+faceNECount;  
             node0 = (periodicLocal[0] && !i0) ? (NE1-1)*N0t :  i0 + periodicLocal[0] + (NE1-1)*N0t;  
   
       out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;  
             out->FaceElements->Tag[k]   = 20;  
             out->FaceElements->Color[k] = i0%2+2;  
             out->FaceElements->Dom[k] = ELEMENT_INTERNAL;  
   
             if (useElementsOnFace) {  
                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;  
                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;  
                 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;  
                 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
             } else {  
                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;  
                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;  
             }  
         }  
       if( boundaryLeft ){  
           out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;  
           if( periodicLocal[0] )  
               out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;  
         }  
       if( boundaryRight ){  
           out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;  
             out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;  
         }  
   
       totalNECount += numElementsLocal;  
       faceNECount  += numElementsLocal;  
     }  
   /* **  elements on boundary 001 (x0=0): */  
   if( domLeft && !periodic[0] )  
     {  
 #pragma omp parallel for private(i1,k,node0)  
         for (i1=0;i1<NE1;i1++) {  
             k=i1+faceNECount;  
             node0=i1*N0t;  
   
       out->FaceElements->Id[k]=faceStart + i1;  
             out->FaceElements->Tag[k]   = 1;  
             out->FaceElements->Color[k] = i1%2+4;  
             out->FaceElements->Dom[k] = ELEMENT_INTERNAL;  
   
             if (useElementsOnFace) {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t+1;  
             } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
   
             }  
         }  
         totalNECount+=NE1;  
         faceNECount+=NE1;  
         }  
         /* **  elements on boundary 002 (x0=1): */  
         if( domRight && !periodic[0])  
         {  
 #pragma omp parallel for private(i1,k,node0)  
         for (i1=0;i1<NE1;i1++) {  
             k=i1+faceNECount;  
             node0=(N0t-2)+i1*N0t;  
   
       out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;  
             out->FaceElements->Tag[k]   = 2;  
             out->FaceElements->Color[k] = i1%2+6;  
             out->FaceElements->Dom[k] = ELEMENT_INTERNAL;  
   
             if (useElementsOnFace) {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
             } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;  
             }  
         }  
         totalNECount+=NE1;  
         faceNECount+=NE1;        
     }  
   out->FaceElements->minColor = 0;  
   out->FaceElements->maxColor = 7;  
   
     out->FaceElements->numElements = faceNECount;  
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
   
   /* setup distribution info for other elements */  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     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;  
   }  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out);  
   if( !Finley_MPI_noError(mpi_info) )  
   {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
   
302    /* free up memory */    /* free up memory */
303    Paso_MPIInfo_dealloc( mpi_info );    Paso_MPIInfo_free( mpi_info );  
   
     //print_mesh_statistics( out, TRUE );  
   
304    #ifdef Finley_TRACE    #ifdef Finley_TRACE
305    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
306    #endif    #endif
   if (Finley_noError()) {  
        if (! Finley_Mesh_isPrepared(out) ) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   }  
   return out;  
 }  
 #endif  
   
   
307    
308      return out;
309    }

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

  ViewVC Help
Powered by ViewVC 1.1.26