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

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

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

revision 751 by bcumming, Mon Jun 26 01:46:34 2006 UTC revision 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC
# Line 43  static index_t domain_MODdim( index_t ra Line 43  static index_t domain_MODdim( index_t ra
43    
44    
45  /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */  /* 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... */  
46  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 )  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 )
47  {  {
48    index_t i0;    index_t i0;
# Line 122  static void domain_calculateDimension( i Line 121  static void domain_calculateDimension( i
121  }  }
122  #endif  #endif
123    
124    #ifdef PASO_MPI
125    Finley_Mesh* Finley_RectangularMesh_Rec4_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
126    #else
127  Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)  Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
128  #ifndef PASO_MPI  #endif
129  {  {
130    dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;    dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;
131    index_t k,node0;    index_t k,node0;
# Line 160  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 162  Finley_Mesh* Finley_RectangularMesh_Rec4
162    /*  allocate mesh: */    /*  allocate mesh: */
163        
164    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
165    #ifdef PASO_MPI
166      out=Finley_Mesh_alloc(name,2,order,mpi_info);
167    
168      if (! Finley_noError()) return NULL;
169    
170      out->Elements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
171      if (useElementsOnFace) {
172         out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order,mpi_info);
173         out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order,mpi_info);
174      } else {
175         out->FaceElements=Finley_ElementFile_alloc(Line2,out->order,mpi_info);
176         out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order,mpi_info);
177      }
178      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
179    
180      if (! Finley_noError()) {
181          Finley_Mesh_dealloc(out);
182          return NULL;
183      }
184      
185      /*  allocate tables: */
186      Finley_NodeFile_allocTable(out->Nodes,N0*N1);
187      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, N0*N1, 0, 0);
188      Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
189      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
190    #else
191    out=Finley_Mesh_alloc(name,2,order);    out=Finley_Mesh_alloc(name,2,order);
192    
193    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
# Line 183  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 211  Finley_Mesh* Finley_RectangularMesh_Rec4
211    Finley_NodeFile_allocTable(out->Nodes,N0*N1);    Finley_NodeFile_allocTable(out->Nodes,N0*N1);
212    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
213    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
214    #endif
215    if (! Finley_noError()) {    if (! Finley_noError()) {
216        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
217        return NULL;        return NULL;
218    }    }
219    /*  set nodes: */    /*  set nodes: */
                                                                                                                                                                                                       
220    #pragma omp parallel for private(i0,i1,k)    #pragma omp parallel for private(i0,i1,k)
221    for (i1=0;i1<N1;i1++) {    for (i1=0;i1<N1;i1++) {
222      for (i0=0;i0<N0;i0++) {      for (i0=0;i0<N0;i0++) {
# Line 199  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 226  Finley_Mesh* Finley_RectangularMesh_Rec4
226        out->Nodes->Id[k]=i0+N0*i1;        out->Nodes->Id[k]=i0+N0*i1;
227        out->Nodes->Tag[k]=0;        out->Nodes->Tag[k]=0;
228        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
229    #ifdef PASO_MPI
230          out->Nodes->Dom[k]=NODE_INTERNAL;
231    #endif
232      }      }
233    }    }
234    /* tags for the faces: */    /* tags for the faces: */
# Line 226  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 256  Finley_Mesh* Finley_RectangularMesh_Rec4
256        out->Elements->Id[k]=k;        out->Elements->Id[k]=k;
257        out->Elements->Tag[k]=0;        out->Elements->Tag[k]=0;
258        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
259    #ifdef PASO_MPI
260          out->Elements->Dom[k]=ELEMENT_INTERNAL;
261    #endif
262    
263        out->Elements->Nodes[INDEX2(0,k,4)]=node0;        out->Elements->Nodes[INDEX2(0,k,4)]=node0;
264        out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;        out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
# Line 257  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 290  Finley_Mesh* Finley_RectangularMesh_Rec4
290            out->FaceElements->Id[k]=i0+totalNECount;            out->FaceElements->Id[k]=i0+totalNECount;
291            out->FaceElements->Tag[k]=10;            out->FaceElements->Tag[k]=10;
292            out->FaceElements->Color[k]=i0%2;            out->FaceElements->Color[k]=i0%2;
293    #ifdef PASO_MPI
294                  out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
295    #endif
296        
297            if (useElementsOnFace) {            if (useElementsOnFace) {
298                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 281  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 317  Finley_Mesh* Finley_RectangularMesh_Rec4
317           out->FaceElements->Id[k]=i0+totalNECount;           out->FaceElements->Id[k]=i0+totalNECount;
318           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
319           out->FaceElements->Color[k]=i0%2+2;           out->FaceElements->Color[k]=i0%2+2;
320    #ifdef PASO_MPI
321                  out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
322    #endif
323        
324           if (useElementsOnFace) {           if (useElementsOnFace) {
325               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0+1;
# Line 305  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 344  Finley_Mesh* Finley_RectangularMesh_Rec4
344            out->FaceElements->Id[k]=i1+totalNECount;            out->FaceElements->Id[k]=i1+totalNECount;
345            out->FaceElements->Tag[k]=1;            out->FaceElements->Tag[k]=1;
346            out->FaceElements->Color[k]=i1%2+4;            out->FaceElements->Color[k]=i1%2+4;
347    #ifdef PASO_MPI
348                  out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
349    #endif
350    
351            if (useElementsOnFace) {            if (useElementsOnFace) {
352                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
# Line 330  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 372  Finley_Mesh* Finley_RectangularMesh_Rec4
372            out->FaceElements->Id[k]=i1+totalNECount;            out->FaceElements->Id[k]=i1+totalNECount;
373            out->FaceElements->Tag[k]=2;            out->FaceElements->Tag[k]=2;
374            out->FaceElements->Color[k]=i1%2+6;            out->FaceElements->Color[k]=i1%2+6;
375    #ifdef PASO_MPI
376                  out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
377    #endif
378        
379            if (useElementsOnFace) {            if (useElementsOnFace) {
380               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;               out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
# Line 346  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 391  Finley_Mesh* Finley_RectangularMesh_Rec4
391    }    }
392    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
393    out->FaceElements->maxColor=7;    out->FaceElements->maxColor=7;
   
     
   /*  face elements done: */  
394        
395    #ifdef PASO_MPI
396        Finley_ElementFile_setDomainFlags( out->Elements );
397        Finley_ElementFile_setDomainFlags( out->FaceElements );
398        Finley_ElementFile_setDomainFlags( out->ContactElements );
399        Finley_ElementFile_setDomainFlags( out->Points );
400    
401        /* reorder the degrees of freedom */
402        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
403    #endif
404        
405    /*   condense the nodes: */    /*   condense the nodes: */
     
406    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
407    
408    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
   
409    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out) ;
410    
411    #ifdef Finley_TRACE    #ifdef Finley_TRACE
# Line 368  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 418  Finley_Mesh* Finley_RectangularMesh_Rec4
418    }    }
419    return out;    return out;
420  }  }
421    #ifdef PASO_MPI
422  /*----------------------------------------------------------------------------  Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
          MPI VERSION  
   
   
   ASSUMPTIONS  
   ===========  
   
   - the domain dimension is large enough in the NE0 direction for each domain  
     to be 2 internal nodes wide in that direction. Face element calculation gets  
     buggered up otherwise.  
   - if dividing into two domains with periodic[0]=TRUE , then the global domain  
     must be at least 4 elements along the NE0 direction, otherwise redundant  
     nodes are generated.  
   
   These problems can be avoided by ensuring that numElements[0]/mpiSize >= 2  
   
   if domains are small enough to cause these problems, you probably don't  
   need to solve them in parallel. If you really do, rotate the domain so that the  
   long axis is aligned with NE0.  
 ------------------------------------------------------------------------------*/  
   
   
 #else  
423  {  {
424    dim_t N0,N1,NE0,NE1,i0,i1, NE0_local, NDOF0,NDOF1, NFaceElements, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, faceNECount, totalNECount,M0,M1,numDOFInternal;    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;
425    index_t innerLoop=0;    index_t NUMNODES,k,firstNode=0, DOFcount=0, node0, node1;
426    index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1;    index_t targetDomain=-1, firstNodeConstruct;
427    index_t targetDomain=-1;    index_t *forwardDOF=NULL, *backwardDOF=NULL;
428    index_t *indexForward=NULL;    bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;  
429    Finley_Mesh* out;    Finley_Mesh* out;
430    char name[50];    char name[50];
431    double time0=Finley_timer();    double time0=Finley_timer();
# Line 414  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 441  Finley_Mesh* Finley_RectangularMesh_Rec4
441    if (! Finley_noError())    if (! Finley_noError())
442          return NULL;          return NULL;
443    
444        /* use the serial code to generate the mesh in the 1-CPU case */
445        if( mpi_info->size==1 ){
446            out = Finley_RectangularMesh_Rec4_singleCPU(numElements,Length,periodic,order,useElementsOnFace,mpi_info);
447            return out;
448        }
449    
450    if( mpi_info->rank==0 )    if( mpi_info->rank==0 )
451      domLeft = TRUE;      domLeft = TRUE;
452    if( mpi_info->rank==mpi_info->size-1 )    if( mpi_info->rank==mpi_info->size-1 )
# Line 446  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 479  Finley_Mesh* Finley_RectangularMesh_Rec4
479    } else {    } else {
480        NDOF1=N1-1;        NDOF1=N1-1;
481    }    }
   
   /*  allocate mesh: */  
482        
483        boundaryLeft = !domLeft || periodicLocal[0];
484        boundaryRight = !domRight || periodicLocal[1];
485        N0t = numNodesLocal + boundaryRight + boundaryLeft;
486        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
487        firstNodeConstruct = firstNode - boundaryLeft;
488        firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
489    
490    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
491    out=Finley_Mesh_alloc( name, 2, order, mpi_info );    out=Finley_Mesh_alloc( name, 2, order, mpi_info );
492    
# Line 470  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 508  Finley_Mesh* Finley_RectangularMesh_Rec4
508        
509    /*  allocate tables: */    /*  allocate tables: */
510    Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );    Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
511    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1, (DOFExternal[0]+DOFExternal[1])*NDOF1, 0 );    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1, 2*NDOF1, 0 );
512    Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);    Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
513    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
514    
# Line 478  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 516  Finley_Mesh* Finley_RectangularMesh_Rec4
516        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
517        return NULL;        return NULL;
518    }    }
519    
520    /*  set nodes: */    /*  set nodes: */
521                                                                                                                                                                                                          #pragma omp parallel for private(i0,i1,k)
   //#pragma omp parallel for private(i0,i1,k)  
   if( mpi_info->size==1 )  
     innerLoop = numNodesLocal;  
   else  
     innerLoop = numNodesLocal-periodicLocal[0];  
     
522    for (k=0,i1=0;i1<N1;i1++) {    for (k=0,i1=0;i1<N1;i1++) {
523      for ( i0=0;i0<innerLoop;i0++,k++) {      for ( i0=0;i0<N0t;i0++,k++) {
524        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
525        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
526        out->Nodes->Id[k]=i0+innerLoop*i1;        out->Nodes->Id[k]=k;
527        out->Nodes->Tag[k]=0;        out->Nodes->Tag[k]=0;
528        out->Nodes->degreeOfFreedom[k]=i0%numDOFLocal + numDOFLocal*(i1%NDOF1);        out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
529                out->Nodes->Dom[k]=NODE_INTERNAL;
530      }      }
531    }    }
532        if( boundaryLeft ){
533            for( i1=0; i1<N1; i1++ ){
534                out->Nodes->Dom[N0t*i1] = NODE_EXTERNAL;
535                out->Nodes->Dom[N0t*i1+1] = NODE_BOUNDARY;
536            }
537        }
538        if( boundaryRight ){
539            for( i1=0; i1<N1; i1++ ){
540                out->Nodes->Dom[N0t*(i1+1)-1] = NODE_EXTERNAL;
541                out->Nodes->Dom[N0t*(i1+1)-2] = NODE_BOUNDARY;
542            }
543        }
544        if( periodicLocal[0] ){
545            for( i1=0; i1<N1; i1++ ){
546                out->Nodes->degreeOfFreedom[i1*N0t+2] = out->Nodes->degreeOfFreedom[i1*N0t+1];
547                out->Nodes->Dom[N0t*i1+2] = NODE_BOUNDARY;
548            }
549        }
550            
551    /* tag 'em */    /* tag 'em */
552    for( i0=0; i0<numNodesLocal; i0++ )    for(i0=0; i0<N0t; i0++ )
553    {    {
554      out->Nodes->Tag[i0] += 10;                        // bottom      out->Nodes->Tag[i0] += 10;                        // bottom
555      out->Nodes->Tag[i0+numNodesLocal*(N1-1)] += 20;   // top      out->Nodes->Tag[i0+N0t*(N1-1)] += 20;   // top
556    }    }
557    if( domLeft && !periodicLocal[0] )    if( domLeft && !periodicLocal[0] )
558      for( i1=0; i1<N1; i1++ )      for( i1=0; i1<N1; i1++ )
559        out->Nodes->Tag[i1*numNodesLocal] += 1;         // left        out->Nodes->Tag[i1*N0t] += 1;         // left
560    
561    if( domRight && !periodicLocal[0] )    if( domRight && !periodicLocal[0] )
562      for( i1=0; i1<N1; i1++ )      for( i1=0; i1<N1; i1++ )
563        out->Nodes->Tag[(i1+1)*numNodesLocal-1] += 2;   // right        out->Nodes->Tag[(i1+1)*N0t-1] += 2;   // right
   
   /* external nodes */  
   k = innerLoop*N1;  
   DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;  
   if( mpi_info->size>1 )  
   {  
     /* add forward communication information for boundary nodes */  
     indexForward = TMPMEMALLOC( NDOF1, index_t );  
     if( domInternal || domRight || periodicLocal[0] )  
     {  
         for( int i=0; i<NDOF1; i++ )  
           indexForward[i] = i*numDOFLocal;  
         targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;  
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
     }  
     if( domInternal || domLeft || periodicLocal[1] )  
     {  
       for( int i=0; i<NDOF1; i++ )  
           indexForward[i] = (i+1)*numDOFLocal - 1;  
       targetDomain = (mpi_info->rank+1) % mpi_info->size;  
       Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
     }  
     TMPMEMFREE( indexForward );  
564    
565      /* LHS */      /* form the boudary communication information */
566      if( periodicLocal[0] )      forwardDOF  = MEMALLOC(NDOF1,index_t);
567      {      backwardDOF = MEMALLOC(NDOF1,index_t);
568        for (i1=0; i1<N1; i1++, k++)      if( !(mpi_info->size==2 && periodicLocal[0])){
569        {          if( boundaryLeft  ) {
570          out->Nodes->Coordinates[INDEX2(0,k,2)]=Length[0];              targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
571          out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];              for( i1=0; i1<NDOF1; i1++ ){
572          out->Nodes->Id[k]=k;                  forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];
573          out->Nodes->Tag[k]=0;                  backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];
574          out->Nodes->degreeOfFreedom[k] = (i1 % NDOF1)*numDOFLocal;              }
575        }              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
576        out->Nodes->Tag[k-N1] = 10;              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
577        out->Nodes->Tag[k-1]  = 20;          }
578        for (i1=0; i1<N1; i1++, k++)          if( boundaryRight ) {
579        {              targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
580          out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(NE0-1)/DBLE(NE0)*Length[0];              for( i1=0; i1<NDOF1; i1++ ){
581          out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];                  forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];
582          out->Nodes->Id[k]=k;                  backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];
583          out->Nodes->Tag[k]=0;              }
584          out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;              Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
585        }              Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
586        out->Nodes->Tag[k-N1] = 10;          }
587        out->Nodes->Tag[k-1]  = 20;      } else{
588        DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;          /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
589            targetDomain = 1;
590        targetDomain = mpi_info->size-1;          
591        Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(numNodesLocal*N1) );          for( i1=0; i1<NDOF1; i1++ ){
592      }                  forwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-2];
593      if( mpi_info->rank>0 )              backwardDOF[i1] = out->Nodes->degreeOfFreedom[(i1+1)*N0t-1];
594      {          }
595        for (i1=0; i1<N1; i1++, k++)          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
596        {          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
597          out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];  
598          out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];          for( i1=0; i1<NDOF1; i1++ ){
599          out->Nodes->Id[k]=k;              forwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t+1];
600          out->Nodes->Tag[k]=0;              backwardDOF[i1] = out->Nodes->degreeOfFreedom[i1*N0t];
601          out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;          }
602        }          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, forwardDOF );
603        out->Nodes->Tag[k-N1] = 10;          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, backwardDOF );
604        out->Nodes->Tag[k-1]  = 20;      }
605        DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;      MEMFREE( forwardDOF );
606        MEMFREE( backwardDOF );
607    
608        targetDomain = mpi_info->rank-1;    /* elements: */
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );  
     }  
       
     /* RHS */  
     if( periodicLocal[1] || (mpi_info->rank < mpi_info->size-1) )  
     {  
       for (i1=0; i1<N1; i1++, k++)  
       {  
         out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];  
         out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];  
         out->Nodes->Id[k]=k;  
         out->Nodes->Tag[k]=0;  
         out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;  
       }  
       out->Nodes->Tag[k-N1] = 10;  
       out->Nodes->Tag[k-1]  = 20;  
       DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;  
         
       targetDomain = (mpi_info->rank+1) % mpi_info->size;  
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );  
     }  
   }  
   out->Nodes->degreeOfFreedomDistribution->numInternal = NDOF1*numDOFInternal;  
   out->Nodes->degreeOfFreedomDistribution->numBoundary = out->Nodes->degreeOfFreedomDistribution->numLocal- out->Nodes->degreeOfFreedomDistribution->numInternal;  
   /*   set the elements: */  
609        
610   /* internal */   #pragma omp parallel for private(i0,i1,k,node0)  
611   //#pragma omp parallel for private(i0,i1,k,node0)     for ( i1=0, k=0; i1<NE1; i1++) {
612   for ( i1=0, k=0; i1<NE1; i1++)      for ( i0=0; i0<numElementsLocal; i0++, k++) {
613   {                node0 = (periodicLocal[0] && !i0) ? i1*N0t :  i1*N0t + i0 + periodicLocal[0];
     for ( i0=0; i0<numElementsInternal; i0++, k++)  
     {  
       node0=i0+i1*innerLoop;  
614    
615        out->Elements->Id[k]=k;        out->Elements->Id[k]=k;
616        out->Elements->Tag[k]=0;        out->Elements->Tag[k]=0;
617        out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
618                out->Elements->Dom[k]=ELEMENT_INTERNAL;
619    
620        out->Elements->Nodes[INDEX2(0,k,4)]=node0;        out->Elements->Nodes[INDEX2(0,k,4)]=node0;
621        out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;        out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
622        out->Elements->Nodes[INDEX2(2,k,4)]=node0+innerLoop+1;        out->Elements->Nodes[INDEX2(2,k,4)]=node0+N0t+1;
623        out->Elements->Nodes[INDEX2(3,k,4)]=node0+innerLoop;        out->Elements->Nodes[INDEX2(3,k,4)]=node0+N0t;
   
     }  
   }  
   /* boundary */  
   if( mpi_info->size>1 )  
   {  
     if( domInternal )  
     {  
       /* left */  
       for( i1=0; i1<NE1; i1++, k++ )  
       {  
         node1=numNodesLocal*N1 + i1;  
         node0=i1*numNodesLocal;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
         out->Elements->Nodes[INDEX2(0,k,4)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,4)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;  
         out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;  
       }  
   
       /* right */  
       for( i1=0; i1<NE1; i1++, k++ )  
       {  
         node0 = (i1+1)*numNodesLocal-1;  
         node1 = (numNodesLocal+1)*N1 + i1;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
         out->Elements->Nodes[INDEX2(0,k,4)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,4)]=node1;  
         out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;  
         out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;  
       }  
624      }      }
     else if( domLeft )  
     {  
       /* left */  
       if( periodicLocal[0] )  
       {  
         node1 = numNodesLocal*N1;  
         node0 = innerLoop*N1;  
         for( i1=0; i1<NE1; i1++, k++, node1++, node0++ )  
         {  
           out->Elements->Id[k]=k;  
           out->Elements->Tag[k]=0;  
           out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
           out->Elements->Nodes[INDEX2(0,k,4)]=node1;  
           out->Elements->Nodes[INDEX2(1,k,4)]=node0;  
           out->Elements->Nodes[INDEX2(2,k,4)]=node0+1;  
           out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;  
         }  
       }  
       /* right */  
       for( i1=0; i1<NE1; i1++, k++ )  
       {  
         node0 = (i1+1)*innerLoop-1;  
         node1 = (numNodesLocal+periodicLocal[0])*N1 + i1;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
         out->Elements->Nodes[INDEX2(0,k,4)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,4)]=node1;  
         out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;  
         out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal-periodicLocal[0];  
       }  
     }  
     else if( domRight )  
     {  
       /* left */  
       for( i1=0; i1<NE1; i1++, k++ )  
       {  
         node1=numNodesLocal*N1 + i1;  
         node0=i1*numNodesLocal;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
         out->Elements->Nodes[INDEX2(0,k,4)]=node1;  
         out->Elements->Nodes[INDEX2(1,k,4)]=node0;  
         out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;  
         out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;  
       }  
       /* right */  
       if( periodicLocal[1] )  
       {  
         for( i1=0; i1<NE1; i1++, k++ )  
         {  
           node0=(i1+1)*numNodesLocal-1;  
           node1 = (numNodesLocal+1)*N1 + i1;  
   
           out->Elements->Id[k]=k;  
           out->Elements->Tag[k]=0;  
           out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);  
   
           out->Elements->Nodes[INDEX2(0,k,4)]=node0;  
           out->Elements->Nodes[INDEX2(1,k,4)]=node1;  
           out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;  
           out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;  
         }  
       }  
     }    
625    }    }
626    out->Elements->minColor=0;    out->Elements->minColor=0;
627    out->Elements->maxColor=0; //COLOR_MOD(0)+3*COLOR_MOD(0);    out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
628        if( boundaryLeft )
629    if( domInternal )          for( i1=0; i1<NE1; i1++ )
630    {              out->Elements->Dom[i1*numElementsLocal]=ELEMENT_BOUNDARY;
631      out->Elements->elementDistribution->numInternal = NE1*(numElementsLocal-2);      if( boundaryRight )
632      out->Elements->elementDistribution->numBoundary = NE1*2;          for( i1=0; i1<NE1; i1++ )
633    }              out->Elements->Dom[(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
634    else if( mpi_info->size==1 )          
635    {      Finley_ElementFile_setDomainFlags( out->Elements );
     out->Elements->elementDistribution->numInternal = NE1*numElementsLocal;  
   }  
   else  
   {  
       out->Elements->elementDistribution->numInternal += NE1*(numElementsLocal-1-periodic[0]);  
       out->Elements->elementDistribution->numBoundary += NE1*(1 + periodic[0]);  
   }  
     
   out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;  
636    
637    /*   face elements: */    /*   face elements: */
     
638    if (useElementsOnFace) {    if (useElementsOnFace) {
639       NUMNODES=4;       NUMNODES=4;
640    } else {    } else {
# Line 760  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 643  Finley_Mesh* Finley_RectangularMesh_Rec4
643    totalNECount=numElementsLocal*NE1;    totalNECount=numElementsLocal*NE1;
644    faceNECount=0;    faceNECount=0;
645        
646    /* do the internal elements first */      if( !periodic[1] ){
647    if (!periodic[1])          
648    {        /* elements on boundary 010 (x1=0): */
649       /* elements on boundary 010 (x1=0): */        #pragma omp parallel for private(i0,k,node0)
650       #pragma omp parallel for private(i0,k,node0)        for (i0=0;i0<numElementsLocal;i0++) {
651       for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {                  k=i0+faceNECount;
652            k=i0+faceNECount;                  node0 = (periodicLocal[0] && !i0) ? 0 :  i0 + periodicLocal[0];
653            node0=i0;  
654                      out->FaceElements->Id[k]=i0+totalNECount;
655            out->FaceElements->Id[k]=i0+totalNECount;                  out->FaceElements->Tag[k]   = 10;
656            out->FaceElements->Tag[k]   = 10;                  out->FaceElements->Color[k] = i0%2;
657            out->FaceElements->Color[k] = 0; // i0%2;                  out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
658    /* WARNING */  
659      /* should be using numDOFLocal instead of numNodesLocal for the case of left hand domain with periodic[0]=true */                  if (useElementsOnFace) {
660            if (useElementsOnFace) {                          out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
661                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                          out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
662                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                          out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t+1;
663                out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal+1;                          out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t;
664                out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;                  } else {
665            } else {                          out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
666                out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                          out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
667                out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                  }
668            }        }
669       }        if( boundaryLeft ){
670       totalNECount += numElementsLocal-numNodesExternal;            out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
671       faceNECount  += numElementsLocal-numNodesExternal;            if( periodicLocal[0] )
672                      out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
673       /* **  elements on boundary 020 (x1=1): */          }
674            if( boundaryRight )
675       #pragma omp parallel for private(i0,k,node0)            out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
      for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {  
          k=i0+faceNECount;  
          node0=i0+(NE1-1)*numNodesLocal;  
     
          out->FaceElements->Id[k]=i0+totalNECount;  
          out->FaceElements->Tag[k]   = 20;  
          out->FaceElements->Color[k] = 0; // i0%2+2;  
     
          if (useElementsOnFace) {  
              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;  
              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+numNodesLocal+1;  
              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;  
          }  
      }  
      totalNECount += numElementsLocal-numNodesExternal;  
      faceNECount  += numElementsLocal-numNodesExternal;  
   }  
   if ( !periodic[0])  
   {  
      /* **  elements on boundary 001 (x0=0): */  
      if( domLeft )  
      {  
        #pragma omp parallel for private(i1,k,node0)  
        for (i1=0;i1<NE1;i1++) {  
             k=i1+faceNECount;  
             node0=i1*numNodesLocal;  
   
             out->FaceElements->Id[k]=i1+totalNECount;  
             out->FaceElements->Tag[k]   = 1;  
             out->FaceElements->Color[k] = 0; //i1%2+4;  
   
             if (useElementsOnFace) {  
                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;  
                 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+numNodesLocal+1;  
             } else {  
                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;  
                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
       
             }  
        }  
        totalNECount+=NE1;  
        faceNECount+=NE1;  
      }  
      /* **  elements on boundary 002 (x0=1): */  
      if( domRight )  
      {  
         #pragma omp parallel for private(i1,k,node0)  
         for (i1=0;i1<NE1;i1++) {  
             k=i1+faceNECount;  
             node0=(numNodesLocal-2)+i1*numNodesLocal;  
   
             out->FaceElements->Id[k]=i1+totalNECount;  
             out->FaceElements->Tag[k]   = 2;  
             out->FaceElements->Color[k] = 0; //i1%2+6;  
         
             if (useElementsOnFace) {  
                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;  
                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;  
                 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;  
                 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+numNodesLocal+1;  
             }  
         }  
         totalNECount+=NE1;  
         faceNECount+=NE1;        
      }  
   }  
   
   if( mpi_info->size>1 )  
     {  
         if( !periodic[1] && (domInternal || domRight) )  
         {  
             // bottom left  
             k     = faceNECount;  
             node0 = numNodesLocal*N1;  
   
             out->FaceElements->Id[k]=totalNECount;  
             out->FaceElements->Tag[k]   = 10;  
             out->FaceElements->Color[k] = 0; //i1%2+6;  
   
             if (useElementsOnFace) {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;  
                     out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal;  
                     out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
             } else {  
                     out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
                     out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;  
             }  
             totalNECount++;  
             faceNECount++;  
676    
677              // top left        totalNECount += numElementsLocal;
678              k     = faceNECount;        faceNECount  += numElementsLocal;
679              node0 = (numNodesLocal+1)*N1 - 2;      
680          /* **  elements on boundary 020 (x1=1): */
681    
682          #pragma omp parallel for private(i0,k,node0)
683          for (i0=0;i0<numElementsLocal;i0++) {
684                k=i0+faceNECount;
685                node0 = (periodicLocal[0] && !i0) ? (NE1-1)*N0t :  i0 + periodicLocal[0] + (NE1-1)*N0t;
686    
687              out->FaceElements->Id[k]=totalNECount;              out->FaceElements->Id[k]=i0+totalNECount;
688              out->FaceElements->Tag[k]   = 20;              out->FaceElements->Tag[k]   = 20;
689              out->FaceElements->Color[k] = 0; //i1%2+6;              out->FaceElements->Color[k] = i0%2+2;
690                out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
691    
692              if (useElementsOnFace) {              if (useElementsOnFace) {
693                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);                  out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;
694                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                  out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;
695                      out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;                  out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
696                      out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=numNodesLocal*(N1-2);                  out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
697              } else {              } else {
698                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);                  out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t+1;
699                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                  out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t;
700              }              }
701              totalNECount++;          }
702              faceNECount++;        if( boundaryLeft ){
703              out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
704              if( periodicLocal[0] )
705                  out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
706          }          }
707          if( !periodic[1] && (domInternal || domLeft) )        if( boundaryRight )
708          {            out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
709              // bottom right  
710              k     = faceNECount;        totalNECount += numElementsLocal;
711              node0 = (numNodesLocal+nodesExternal[0])*N1;        faceNECount  += numElementsLocal;
712        }
713              out->FaceElements->Id[k]=totalNECount;    /* **  elements on boundary 001 (x0=0): */
714              out->FaceElements->Tag[k]   = 10;    if( domLeft && !periodic[0] )
715              out->FaceElements->Color[k] = 0; //i1%2+6;      {
716    #pragma omp parallel for private(i1,k,node0)
717            for (i1=0;i1<NE1;i1++) {
718                k=i1+faceNECount;
719                node0=i1*N0t;
720    
721                out->FaceElements->Id[k]=i1+totalNECount;
722                out->FaceElements->Tag[k]   = 1;
723                out->FaceElements->Color[k] = i1%2+4;
724                out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
725    
726              if (useElementsOnFace) {              if (useElementsOnFace) {
727                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
728                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
729                      out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;                      out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
730                      out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=2*numNodesLocal-1;                      out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0t+1;
731              } else {              } else {
732                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0t;
733                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
             }  
             totalNECount++;  
             faceNECount++;  
734    
735              // top right              }
736              k     = faceNECount;          }
737              node0 = (numNodesLocal+1+nodesExternal[0])*N1-2;          totalNECount+=NE1;
738            faceNECount+=NE1;
739              out->FaceElements->Id[k]=totalNECount;          }
740              out->FaceElements->Tag[k]   = 20;          /* **  elements on boundary 002 (x0=1): */
741              out->FaceElements->Color[k] = 0; //i1%2+6;          if( domRight && !periodic[0])
742            {
743    #pragma omp parallel for private(i1,k,node0)
744            for (i1=0;i1<NE1;i1++) {
745                k=i1+faceNECount;
746                node0=(N0t-2)+i1*N0t;
747    
748                out->FaceElements->Id[k]=i1+totalNECount;
749                out->FaceElements->Tag[k]   = 2;
750                out->FaceElements->Color[k] = i1%2+6;
751                out->FaceElements->Dom[k] = ELEMENT_INTERNAL;
752    
753              if (useElementsOnFace) {              if (useElementsOnFace) {
754                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
755                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
756                      out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal*(N1-1)-1;                      out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0t;
757                      out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;                      out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
758              } else {              } else {
759                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                      out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
760                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;                      out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0t+1;
761              }              }
             totalNECount++;  
             faceNECount++;  
762          }          }
763            totalNECount+=NE1;
764            faceNECount+=NE1;      
765      }      }
766    out->FaceElements->minColor = 0;    out->FaceElements->minColor = 0;
767    out->FaceElements->maxColor = 0; //7;    out->FaceElements->maxColor = 7;
768    
769      out->FaceElements->elementDistribution->numInternal = 0;      out->FaceElements->numElements = faceNECount;
770      out->FaceElements->elementDistribution->numBoundary = 0;      Finley_ElementFile_setDomainFlags( out->FaceElements );
   if( domInternal )  
   {  
     if( !periodic[1] )  
     {  
       out->FaceElements->elementDistribution->numInternal = 2*(numElementsLocal-2);  
       out->FaceElements->elementDistribution->numBoundary = 4;  
     }  
   }  
   else if( mpi_info->size==1 )  
   {  
     if( !periodic[0] )  
       out->FaceElements->elementDistribution->numInternal += NE1*2;  
     if( !periodic[1] )  
       out->FaceElements->elementDistribution->numInternal += numElementsLocal*2;  
   }  
   else  
   {  
     if( !periodic[0] )  
       out->FaceElements->elementDistribution->numInternal += NE1;  
     if( !periodic[1] )  
     {  
       out->FaceElements->elementDistribution->numInternal += 2*(numElementsLocal-1-periodic[0]);  
       out->FaceElements->elementDistribution->numBoundary += 2*(1 + periodic[0]);  
     }  
   }  
     
   out->FaceElements->elementDistribution->numLocal = faceNECount;  
771    
   /*  face elements done: */  
       
772    /* setup distribution info for other elements */    /* setup distribution info for other elements */
773    out->ContactElements->elementDistribution->numLocal = out->ContactElements->elementDistribution->numInternal = out->ContactElements->elementDistribution->numInternal = 0;      Finley_ElementFile_setDomainFlags( out->ContactElements );
774    out->Points->elementDistribution->numLocal = out->Points->elementDistribution->numInternal = out->Points->elementDistribution->numInternal = 0;      Finley_ElementFile_setDomainFlags( out->Points );
775      
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds( out );  
776    
777    /* setup the CommBuffer */      /* reorder the degrees of freedom */
778    Finley_NodeDistribution_formCommBuffer( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );      Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
779    if ( !Finley_MPI_noError( mpi_info )) {      
780      if( Finley_noError() )    /*   condense the nodes: */
781        Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );    Finley_Mesh_resolveNodeIds(out);
782      if( !Finley_MPI_noError(mpi_info) )
783      {
784      Paso_MPIInfo_dealloc( mpi_info );      Paso_MPIInfo_dealloc( mpi_info );
785      Finley_Mesh_dealloc(out);      Finley_Mesh_dealloc(out);
786      return NULL;      return NULL;
787    }    }
   
   Finley_NodeDistribution_calculateIndexExternal( out->Nodes->degreeOfFreedomDistribution, out->Nodes->CommBuffer );  
788    
789    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
790    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out);
   
   #ifdef Finley_TRACE  
   printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);  
   #endif  
   
791    if( !Finley_MPI_noError(mpi_info) )    if( !Finley_MPI_noError(mpi_info) )
792    {    {
     if( Finley_noError() )  
       Finley_setError( PASO_MPI_ERROR, "Error on another MPI process" );  
793      Paso_MPIInfo_dealloc( mpi_info );      Paso_MPIInfo_dealloc( mpi_info );
794      Finley_Mesh_dealloc(out);      Finley_Mesh_dealloc(out);
795      return NULL;      return NULL;
796    }    }
797      
798      /* free up memory */    /* free up memory */
799    Paso_MPIInfo_dealloc( mpi_info );    Paso_MPIInfo_dealloc( mpi_info );
800    
801    return out;      //print_mesh_statistics( out, TRUE );
802    
803      #ifdef Finley_TRACE
804      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
805      #endif
806    
807      return out;
808  }  }
809  #endif  #endif
810    

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

  ViewVC Help
Powered by ViewVC 1.1.26