/[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 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 730 by bcumming, Mon May 15 04:03:49 2006 UTC
# Line 29  Line 29 
29    
30  /**************************************************************/  /**************************************************************/
31    
32  Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {  #ifdef PASO_MPI
33    static void Finley_Mesh_printDistributed( Finley_Mesh *in )
34    {
35      dim_t k, i0, NUMNODES;
36      Finley_Mesh *out=in;
37    
38      NUMNODES = in->FaceElements->ReferenceElement->Type->numNodes;
39    
40      printf( "\n============NODES=============\n" );
41      for( k=0; k<in->Nodes->numNodes; k++ )
42        printf( "\tId %d\tDOF %d\tcoord [%3g%3g]\n", out->Nodes->Id[k], out->Nodes->degreeOfFreedom[k] , out->Nodes->Coordinates[INDEX2(0,k,2)], out->Nodes->Coordinates[INDEX2(1,k,2)] );
43        for( k=0; k<out->Nodes->degreeOfFreedomDistribution->numNeighbours; k++ )
44        {
45          if( out->Nodes->degreeOfFreedomDistribution->neighbours[k]>=0 )
46          {
47            printf( "\t%d boundary DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward );
48            for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numForward; i0++ )
49              printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexForward[i0] );
50            printf("} to %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
51            printf( "\t%d external DOF { ", out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward );
52            for( i0=0; i0<out->Nodes->degreeOfFreedomDistribution->edges[k]->numBackward; i0++ )
53              printf( "%d ", out->Nodes->degreeOfFreedomDistribution->edges[k]->indexBackward[i0] );
54            printf("} from %d\n", out->Nodes->degreeOfFreedomDistribution->neighbours[k] );
55          }
56        }
57    
58      printf( "\n============ELEMENTS=============\n");
59      for( k=0; k<in->Elements->elementDistribution->numInternal; k++ )
60      {
61        printf( "I\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,4)], out->Elements->Nodes[INDEX2(1,k,4)], out->Elements->Nodes[INDEX2(2,k,4)], out->Elements->Nodes[INDEX2(3,k,4)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(3,k,4)]] );
62      }
63    
64      for( k=in->Elements->elementDistribution->numInternal; k<in->Elements->elementDistribution->numInternal+in->Elements->elementDistribution->numBoundary; k++ )
65      {
66        printf( "B\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->Elements->Id[k], out->Elements->Nodes[INDEX2(0,k,4)], out->Elements->Nodes[INDEX2(1,k,4)], out->Elements->Nodes[INDEX2(2,k,4)], out->Elements->Nodes[INDEX2(3,k,4)], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(0,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(1,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(2,k,4)]], out->Nodes->degreeOfFreedom[out->Elements->Nodes[INDEX2(3,k,4)]] );
67      }
68    
69      for( k=in->FaceElements->elementDistribution->numInternal; k<in->FaceElements->elementDistribution->numInternal+in->FaceElements->elementDistribution->numBoundary; k++ )
70      {
71        if( NUMNODES==4 )
72          printf( "F\tId %d : nodes [%d %d %d %d]->DOF [%d %d %d %d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]] );
73        else
74          printf( "F\tId %d : nodes [%d %d]->DOF [%d %d]\n", out->FaceElements->Id[k], out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)], out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]], out->Nodes->degreeOfFreedom[out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]] );
75      }
76    
77    }
78    #endif
79    
80    
81    
82    /* get the number of nodes/elements for domain with rank=rank, of size processors
83       where n is the total number of nodes/elements in the global domain */
84    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
85    {
86      rank = size-rank-1;
87    
88      if( rank < n%size )
89        return (index_t)floor(n/size)+1;
90      return (index_t)floor(n/size);
91    }
92    
93    
94    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
95    /* A bit messy, but it only has to be done once... */
96    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 )
97    {
98      index_t i0;
99      dim_t numNodesGlobal = numElementsGlobal+1;
100    
101      (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
102      
103      numElementsLocal[0] = numNodesLocal[0]+1;
104      periodicLocal[0] = periodicLocal[1] = FALSE;
105      nodesExternal[0] = nodesExternal[1] = 1;
106      if( periodic )
107      {
108        if( size==1 )
109        {
110          numElementsLocal[0] = numElementsGlobal;
111          nodesExternal[0] = nodesExternal[1] = 0;
112          periodicLocal[0] = periodicLocal[1] = TRUE;
113        }
114        else
115        {
116          if( rank==0 )
117          {
118            periodicLocal[0] = TRUE;
119            numNodesLocal[0]++;
120          }
121          if( rank==(size-1) )
122          {      
123            periodicLocal[1] = TRUE;
124            numNodesLocal[0]--;
125            numElementsLocal[0]--;
126          }
127        }
128      }
129      else if( !periodic )
130      {
131        if( rank==0 ){
132          nodesExternal[0]--;
133          numElementsLocal[0]--;
134        }
135        if( rank==(size-1) )
136        {
137          nodesExternal[1]--;
138          numElementsLocal[0]--;
139        }
140      }
141      numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
142      
143      numElementsInternal[0] = numElementsLocal[0];
144      if( (rank==0) && (rank==size-1) );
145        
146      else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
147          numElementsInternal[0] -= 1;
148      else
149        numElementsInternal[0] -= 2;
150    
151      firstNode[0] = 0;
152      for( i0=0; i0<rank; i0++ )
153        firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
154    
155      numDOFLocal[0] = numNodesLocal[0];
156      if( periodicLocal[0] )
157      {
158        numDOFLocal[0]--;
159      }
160      DOFExternal[0] = nodesExternal[0];
161      DOFExternal[1] = nodesExternal[1];
162    
163      /* some debugging printf statements */
164      //printf( "rank/size = %d/%d\nNodes : %d Local, %d External[%d %d], First = %d\nElements : %d Local\nDOF : %d Local, External [%d %d]\nperiodicLocal [%d %d]\n\n", rank, size, *numNodesLocal, *numNodesExternal, nodesExternal[0], nodesExternal[1], *firstNode, *numElementsLocal, *numDOFLocal, DOFExternal[0], DOFExternal[1], periodicLocal[0], periodicLocal[1] );
165    }
166    
167    
168    Finley_Mesh* Finley_RectangularMesh_Rec4(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
169    #ifndef PASO_MPI
170    {
171    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;
172    index_t k,node0;    index_t k,node0;
173    Finley_Mesh* out;    Finley_Mesh* out;
# Line 66  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 204  Finley_Mesh* Finley_RectangularMesh_Rec4
204        
205    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
206    out=Finley_Mesh_alloc(name,2,order);    out=Finley_Mesh_alloc(name,2,order);
207    
208    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
209    
210    out->Elements=Finley_ElementFile_alloc(Rec4,out->order);    out->Elements=Finley_ElementFile_alloc(Rec4,out->order);
# Line 77  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 216  Finley_Mesh* Finley_RectangularMesh_Rec4
216       out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);       out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order);
217    }    }
218    out->Points=Finley_ElementFile_alloc(Point1,out->order);    out->Points=Finley_ElementFile_alloc(Point1,out->order);
219    
220    if (! Finley_noError()) {    if (! Finley_noError()) {
221        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
222        return NULL;        return NULL;
223    }    }
224        
225    /*  allocate tables: */    /*  allocate tables: */
     
226    Finley_NodeFile_allocTable(out->Nodes,N0*N1);    Finley_NodeFile_allocTable(out->Nodes,N0*N1);
227    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
228    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
229    
230    if (! Finley_noError()) {    if (! Finley_noError()) {
231        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
232        return NULL;        return NULL;
# Line 250  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 390  Finley_Mesh* Finley_RectangularMesh_Rec4
390    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
391    out->FaceElements->maxColor=7;    out->FaceElements->maxColor=7;
392    
393      
394    /*  face elements done: */    /*  face elements done: */
395        
396    /*   condense the nodes: */    /*   condense the nodes: */
# Line 270  Finley_Mesh* Finley_RectangularMesh_Rec4 Line 411  Finley_Mesh* Finley_RectangularMesh_Rec4
411    }    }
412    return out;    return out;
413  }  }
414    
415    /*----------------------------------------------------------------------------
416             MPI VERSION
417    
418    
419      ASSUMPTIONS
420      ===========
421    
422      - the domain dimension is large enough in the NE0 direction for each domain
423        to be 2 internal nodes wide in that direction. Face element calculation gets
424        buggered up otherwise.
425      - if dividing into two domains with periodic[0]=TRUE , then the global domain
426        must be at least 4 elements along the NE0 direction, otherwise redundant
427        nodes are generated.
428    
429      These problems can be avoided by ensuring that numElements[0]/mpiSize >= 2
430    
431      if domains are small enough to cause these problems, you probably don't
432      need to solve them in parallel. If you really do, rotate the domain so that the
433      long axis is aligned with NE0.
434    ------------------------------------------------------------------------------*/
435    
436    
437    #else
438    {
439      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;
440      index_t innerLoop=0;
441      index_t NUMNODES,k,firstNode=0, DOFcount=0, forwardDOF[2], backwardDOF[2], node0, node1;
442      index_t targetDomain=-1;
443      index_t *indexForward=NULL;
444      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE;
445      Finley_Mesh* out;
446      char name[50];
447      double time0=Finley_timer();
448    
449      NE0=MAX(1,numElements[0]);
450      NE1=MAX(1,numElements[1]);
451      N0=NE0+1;
452      N1=NE1+1;
453      Paso_MPIInfo *mpi_info = NULL;
454    
455      /* get MPI information */
456      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
457      if (! Finley_noError())
458            return NULL;
459    
460      if( mpi_info->rank==0 )
461        domLeft = TRUE;
462      if( mpi_info->rank==mpi_info->size-1 )
463        domRight = TRUE;
464      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
465        domInternal = TRUE;
466    
467      /* dimensions of the local subdomain */
468      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );  
469    
470      if (N0<=N1) {
471         M0=1;
472         M1=N0;
473      } else {
474         M0=N1;
475         M1=1;
476      }
477    
478      NFaceElements=0;
479      if (!periodic[0])
480      {
481          if( domLeft )
482            NFaceElements+=NE1;
483          if( domRight )
484            NFaceElements+=NE1;
485      }
486      if (!periodic[1]) {
487          NDOF1=N1;
488          NFaceElements+=2*NE0;
489      } else {
490          NDOF1=N1-1;
491      }
492    
493      /*  allocate mesh: */
494      
495      sprintf(name,"Rectangular %d x %d mesh",N0,N1);
496      out=Finley_Mesh_alloc( name, 2, order, mpi_info );
497    
498      if (! Finley_noError()) return NULL;
499    
500      out->Elements=Finley_ElementFile_alloc( Rec4, out->order, mpi_info );
501      if (useElementsOnFace) {
502         out->FaceElements=Finley_ElementFile_alloc(Rec4Face,out->order, mpi_info );
503         out->ContactElements=Finley_ElementFile_alloc(Rec4Face_Contact,out->order, mpi_info );
504      } else {
505         out->FaceElements=Finley_ElementFile_alloc(Line2,out->order, mpi_info );
506         out->ContactElements=Finley_ElementFile_alloc(Line2_Contact,out->order, mpi_info );
507      }
508      out->Points=Finley_ElementFile_alloc(Point1,out->order, mpi_info );
509      if (! Finley_noError()) {
510          Finley_Mesh_dealloc(out);
511          return NULL;
512      }
513      
514      /*  allocate tables: */
515      Finley_NodeFile_allocTable( out->Nodes, (numNodesLocal + numNodesExternal)*N1 );
516      Finley_NodeDistibution_allocTable( out->Nodes->degreeOfFreedomDistribution, numNodesLocal*NDOF1, numNodesExternal*NDOF1, 0 );
517      Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
518      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
519    
520      if (! Finley_noError()) {
521          Finley_Mesh_dealloc(out);
522          return NULL;
523      }
524      /*  set nodes: */
525                                                                                                                                                                                                        
526      //#pragma omp parallel for private(i0,i1,k)
527      if( mpi_info->size==1 )
528        innerLoop = numNodesLocal;
529      else
530        innerLoop = numNodesLocal-periodicLocal[0];
531      
532      for (k=0,i1=0;i1<N1;i1++) {
533        for ( i0=0;i0<innerLoop;i0++,k++) {
534          out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0+firstNode)/DBLE(N0-1)*Length[0];
535          out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
536          out->Nodes->Id[k]=i0+innerLoop*i1;
537          out->Nodes->Tag[k]=0;
538          out->Nodes->degreeOfFreedom[k]=i0%numDOFLocal + numDOFLocal*(i1%NDOF1);
539        }
540      }
541    
542      /* tag 'em */
543      for( i0=0; i0<numNodesLocal; i0++ )
544      {
545        out->Nodes->Tag[i0] += 10;                        // bottom
546        out->Nodes->Tag[i0+numNodesLocal*(N1-1)] += 20;   // top
547      }
548      if( domLeft && !periodicLocal[0] )
549        for( i1=0; i1<N1; i1++ )
550          out->Nodes->Tag[i1*numNodesLocal] += 1;         // left
551    
552      if( domRight && !periodicLocal[0] )
553        for( i1=0; i1<N1; i1++ )
554          out->Nodes->Tag[(i1+1)*numNodesLocal-1] += 2;   // right
555    
556      /* external nodes */
557      k = innerLoop*N1;
558      DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
559      if( mpi_info->size>1 )
560      {
561        /* add forward communication information for boundary nodes */
562        indexForward = TMPMEMALLOC( NDOF1, index_t );
563        if( domInternal || domRight || periodicLocal[0] )
564        {
565            for( int i=0; i<NDOF1; i++ )
566              indexForward[i] = i*numDOFLocal;
567            targetDomain = mpi_info->rank-1>=0 ? mpi_info->rank-1 : mpi_info->size-1;
568            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
569        }
570        if( domInternal || domLeft || periodicLocal[1] )
571        {
572          for( int i=0; i<NDOF1; i++ )
573              indexForward[i] = (i+1)*numDOFLocal - 1;
574          targetDomain = (mpi_info->rank+1) % mpi_info->size;
575          Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
576        }
577        TMPMEMFREE( indexForward );
578    
579        /* LHS */
580        if( periodicLocal[0] )
581        {
582          for (i1=0; i1<N1; i1++, k++)
583          {
584            out->Nodes->Coordinates[INDEX2(0,k,2)]=Length[0];
585            out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
586            out->Nodes->Id[k]=k;
587            out->Nodes->Tag[k]=0;
588            out->Nodes->degreeOfFreedom[k] = (i1 % NDOF1)*numDOFLocal;
589          }
590          out->Nodes->Tag[k-N1] = 10;
591          out->Nodes->Tag[k-1]  = 20;
592          for (i1=0; i1<N1; i1++, k++)
593          {
594            out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(NE0-1)/DBLE(NE0)*Length[0];
595            out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
596            out->Nodes->Id[k]=k;
597            out->Nodes->Tag[k]=0;
598            out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
599          }
600          out->Nodes->Tag[k-N1] = 10;
601          out->Nodes->Tag[k-1]  = 20;
602          DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
603    
604          targetDomain = mpi_info->size-1;
605          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(numNodesLocal*N1) );
606        }    
607        if( mpi_info->rank>0 )
608        {
609          for (i1=0; i1<N1; i1++, k++)
610          {
611            out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode-1)/DBLE(N0-1)*Length[0];
612            out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
613            out->Nodes->Id[k]=k;
614            out->Nodes->Tag[k]=0;
615            out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
616          }
617          out->Nodes->Tag[k-N1] = 10;
618          out->Nodes->Tag[k-1]  = 20;
619          DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
620    
621          targetDomain = mpi_info->rank-1;
622          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
623        }
624        
625        /* RHS */
626        if( periodicLocal[1] || (mpi_info->rank < mpi_info->size-1) )
627        {
628          for (i1=0; i1<N1; i1++, k++)
629          {
630            out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(firstNode+numNodesLocal-periodicLocal[0])/DBLE(N0-1)*Length[0];
631            out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
632            out->Nodes->Id[k]=k;
633            out->Nodes->Tag[k]=0;
634            out->Nodes->degreeOfFreedom[k] = DOFcount + i1%NDOF1;
635          }
636          out->Nodes->Tag[k-N1] = 10;
637          out->Nodes->Tag[k-1]  = 20;
638          DOFcount = out->Nodes->degreeOfFreedom[k-1]+1;
639          
640          targetDomain = (mpi_info->rank+1) % mpi_info->size;
641          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, out->Nodes->degreeOfFreedom+(k-N1) );
642        }
643      }
644    
645      /*   set the elements: */
646      
647     /* internal */
648     //#pragma omp parallel for private(i0,i1,k,node0)  
649     for ( i1=0, k=0; i1<NE1; i1++)
650     {  
651        for ( i0=0; i0<numElementsInternal; i0++, k++)
652        {
653          node0=i0+i1*innerLoop;
654    
655          out->Elements->Id[k]=k;
656          out->Elements->Tag[k]=0;
657          out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
658    
659          out->Elements->Nodes[INDEX2(0,k,4)]=node0;
660          out->Elements->Nodes[INDEX2(1,k,4)]=node0+1;
661          out->Elements->Nodes[INDEX2(2,k,4)]=node0+innerLoop+1;
662          out->Elements->Nodes[INDEX2(3,k,4)]=node0+innerLoop;
663    
664        }
665      }
666      /* boundary */
667      if( mpi_info->size>1 )
668      {
669        if( domInternal )
670        {
671          /* left */
672          for( i1=0; i1<NE1; i1++, k++ )
673          {
674            node1=numNodesLocal*N1 + i1;
675            node0=i1*numNodesLocal;
676    
677            out->Elements->Id[k]=k;
678            out->Elements->Tag[k]=0;
679            out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
680    
681            out->Elements->Nodes[INDEX2(0,k,4)]=node1;
682            out->Elements->Nodes[INDEX2(1,k,4)]=node0;
683            out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
684            out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
685          }
686    
687          /* right */
688          for( i1=0; i1<NE1; i1++, k++ )
689          {
690            node0 = (i1+1)*numNodesLocal-1;
691            node1 = (numNodesLocal+1)*N1 + i1;
692    
693            out->Elements->Id[k]=k;
694            out->Elements->Tag[k]=0;
695            out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
696    
697            out->Elements->Nodes[INDEX2(0,k,4)]=node0;
698            out->Elements->Nodes[INDEX2(1,k,4)]=node1;
699            out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
700            out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
701          }
702        }
703        else if( domLeft )
704        {
705          /* left */
706          if( periodicLocal[0] )
707          {
708            node1 = numNodesLocal*N1;
709            node0 = innerLoop*N1;
710            for( i1=0; i1<NE1; i1++, k++, node1++, node0++ )
711            {
712              out->Elements->Id[k]=k;
713              out->Elements->Tag[k]=0;
714              out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
715    
716              out->Elements->Nodes[INDEX2(0,k,4)]=node1;
717              out->Elements->Nodes[INDEX2(1,k,4)]=node0;
718              out->Elements->Nodes[INDEX2(2,k,4)]=node0+1;
719              out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
720            }
721          }
722          /* right */
723          for( i1=0; i1<NE1; i1++, k++ )
724          {
725            node0 = (i1+1)*innerLoop-1;
726            node1 = (numNodesLocal+periodicLocal[0])*N1 + i1;
727    
728            out->Elements->Id[k]=k;
729            out->Elements->Tag[k]=0;
730            out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
731    
732            out->Elements->Nodes[INDEX2(0,k,4)]=node0;
733            out->Elements->Nodes[INDEX2(1,k,4)]=node1;
734            out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
735            out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal-periodicLocal[0];
736          }
737        }
738        else if( domRight )
739        {
740          /* left */
741          for( i1=0; i1<NE1; i1++, k++ )
742          {
743            node1=numNodesLocal*N1 + i1;
744            node0=i1*numNodesLocal;
745    
746            out->Elements->Id[k]=k;
747            out->Elements->Tag[k]=0;
748            out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
749    
750            out->Elements->Nodes[INDEX2(0,k,4)]=node1;
751            out->Elements->Nodes[INDEX2(1,k,4)]=node0;
752            out->Elements->Nodes[INDEX2(2,k,4)]=node0+numNodesLocal;
753            out->Elements->Nodes[INDEX2(3,k,4)]=node1+1;
754          }
755          /* right */
756          if( periodicLocal[1] )
757          {
758            for( i1=0; i1<NE1; i1++, k++ )
759            {
760              node0=(i1+1)*numNodesLocal-1;
761              node1 = (numNodesLocal+1)*N1 + i1;
762    
763              out->Elements->Id[k]=k;
764              out->Elements->Tag[k]=0;
765              out->Elements->Color[k]=0; //COLOR_MOD(i0)+3*COLOR_MOD(i1);
766    
767              out->Elements->Nodes[INDEX2(0,k,4)]=node0;
768              out->Elements->Nodes[INDEX2(1,k,4)]=node1;
769              out->Elements->Nodes[INDEX2(2,k,4)]=node1+1;
770              out->Elements->Nodes[INDEX2(3,k,4)]=node0+numNodesLocal;
771            }
772          }
773        }  
774      }
775      out->Elements->minColor=0;
776      out->Elements->maxColor=0; //COLOR_MOD(0)+3*COLOR_MOD(0);
777    
778      if( domInternal )
779      {
780        out->Elements->elementDistribution->numInternal = NE1*(numElementsLocal-2);
781        out->Elements->elementDistribution->numBoundary = NE1*2;
782      }
783      else if( mpi_info->size==1 )
784      {
785        out->Elements->elementDistribution->numInternal = NE1*numElementsLocal;
786      }
787      else
788      {
789          out->Elements->elementDistribution->numInternal += NE1*(numElementsLocal-1-periodic[0]);
790          out->Elements->elementDistribution->numBoundary += NE1*(1 + periodic[0]);
791      }
792      
793      out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
794    
795      /*   face elements: */
796      
797      if (useElementsOnFace) {
798         NUMNODES=4;
799      } else {
800         NUMNODES=2;
801      }
802      totalNECount=numElementsLocal*NE1;
803      faceNECount=0;
804      
805      /* do the internal elements first */
806      if (!periodic[1])
807      {
808         /* elements on boundary 010 (x1=0): */
809         #pragma omp parallel for private(i0,k,node0)
810         for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
811              k=i0+faceNECount;
812              node0=i0;
813      
814              out->FaceElements->Id[k]=i0+totalNECount;
815              out->FaceElements->Tag[k]   = 10;
816              out->FaceElements->Color[k] = 0; // i0%2;
817      
818              if (useElementsOnFace) {
819                  out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
820                  out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
821                  out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal+1;
822                  out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal;
823              } else {
824                  out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
825                  out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
826              }
827         }
828         totalNECount += numElementsLocal-numNodesExternal;
829         faceNECount  += numElementsLocal-numNodesExternal;
830        
831         /* **  elements on boundary 020 (x1=1): */
832      
833         #pragma omp parallel for private(i0,k,node0)
834         for (i0=0;i0<numElementsLocal-numNodesExternal;i0++) {
835             k=i0+faceNECount;
836             node0=i0+(NE1-1)*numNodesLocal;
837      
838             out->FaceElements->Id[k]=i0+totalNECount;
839             out->FaceElements->Tag[k]   = 20;
840             out->FaceElements->Color[k] = 0; // i0%2+2;
841      
842             if (useElementsOnFace) {
843                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
844                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
845                 out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
846                 out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
847             } else {
848                 out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal+1;
849                 out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal;
850             }
851         }
852         totalNECount += numElementsLocal-numNodesExternal;
853         faceNECount  += numElementsLocal-numNodesExternal;
854      }
855      if ( !periodic[0])
856      {
857         /* **  elements on boundary 001 (x0=0): */
858         if( domLeft )
859         {
860           #pragma omp parallel for private(i1,k,node0)
861           for (i1=0;i1<NE1;i1++) {
862                k=i1+faceNECount;
863                node0=i1*numNodesLocal;
864    
865                out->FaceElements->Id[k]=i1+totalNECount;
866                out->FaceElements->Tag[k]   = 1;
867                out->FaceElements->Color[k] = 0; //i1%2+4;
868    
869                if (useElementsOnFace) {
870                    out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
871                    out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
872                    out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
873                    out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+numNodesLocal+1;
874                } else {
875                    out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+numNodesLocal;
876                    out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
877        
878                }
879           }
880           totalNECount+=NE1;
881           faceNECount+=NE1;
882         }
883         /* **  elements on boundary 002 (x0=1): */
884         if( domRight )
885         {
886            #pragma omp parallel for private(i1,k,node0)
887            for (i1=0;i1<NE1;i1++) {
888                k=i1+faceNECount;
889                node0=(numNodesLocal-2)+i1*numNodesLocal;
890    
891                out->FaceElements->Id[k]=i1+totalNECount;
892                out->FaceElements->Tag[k]   = 2;
893                out->FaceElements->Color[k] = 0; //i1%2+6;
894          
895                if (useElementsOnFace) {
896                    out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
897                    out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
898                    out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+numNodesLocal;
899                    out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
900                } else {
901                    out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
902                    out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+numNodesLocal+1;
903                }
904            }
905            totalNECount+=NE1;
906            faceNECount+=NE1;      
907         }
908      }
909    
910      if( mpi_info->size>1 )
911    {
912      if( !periodic[1] && (domInternal || domRight) )
913      {
914        // bottom left
915        k     = faceNECount;
916        node0 = numNodesLocal*N1;
917    
918        out->FaceElements->Id[k]=totalNECount;
919        out->FaceElements->Tag[k]   = 10;
920        out->FaceElements->Color[k] = 0; //i1%2+6;
921    
922        if (useElementsOnFace) {
923            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
924            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
925            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal;
926            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;
927        } else {
928            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
929            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=0;
930        }
931        totalNECount++;
932        faceNECount++;
933    
934        // top left
935        k     = faceNECount;
936        node0 = (numNodesLocal+1)*N1 - 2;
937    
938        out->FaceElements->Id[k]=totalNECount;
939        out->FaceElements->Tag[k]   = 20;
940        out->FaceElements->Color[k] = 0; //i1%2+6;
941    
942        if (useElementsOnFace) {
943            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
944            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
945            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;
946            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=numNodesLocal*(N1-2);
947        } else {
948            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal*(N1-1);
949            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;
950        }
951        totalNECount++;
952        faceNECount++;
953      }
954      if( !periodic[1] && (domInternal || domLeft) )
955      {
956        // bottom right
957        k     = faceNECount;
958        node0 = (numNodesLocal+nodesExternal[0])*N1;
959    
960        out->FaceElements->Id[k]=totalNECount;
961        out->FaceElements->Tag[k]   = 10;
962        out->FaceElements->Color[k] = 0; //i1%2+6;
963    
964        if (useElementsOnFace) {
965            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
966            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
967            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;
968            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=2*numNodesLocal-1;
969        } else {
970            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=numNodesLocal-1;
971            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;
972        }
973        totalNECount++;
974        faceNECount++;
975    
976        // top right
977        k     = faceNECount;
978        node0 = (numNodesLocal+1+nodesExternal[0])*N1-2;
979    
980        out->FaceElements->Id[k]=totalNECount;
981        out->FaceElements->Tag[k]   = 20;
982        out->FaceElements->Color[k] = 0; //i1%2+6;
983    
984        if (useElementsOnFace) {
985            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
986            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
987            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=numNodesLocal*(N1-1)-1;
988            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;
989        } else {
990            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
991            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=numNodesLocal*N1-1;
992        }
993        totalNECount++;
994        faceNECount++;
995      }
996    }
997      out->FaceElements->minColor = 0;
998      out->FaceElements->maxColor = 0; //7;
999    
1000      if( domInternal )
1001      {
1002        if( !periodic[1] )
1003        {
1004          out->FaceElements->elementDistribution->numInternal = 2*(numElementsLocal-2);
1005          out->FaceElements->elementDistribution->numBoundary = 4;
1006        }
1007      }
1008      else if( mpi_info->size==1 )
1009      {
1010        if( !periodic[0] )
1011          out->FaceElements->elementDistribution->numInternal += NE1*2;
1012        if( !periodic[1] )
1013          out->FaceElements->elementDistribution->numInternal += numElementsLocal*2;
1014      }
1015      else
1016      {
1017        if( !periodic[0] )
1018          out->FaceElements->elementDistribution->numInternal += NE1;
1019        if( !periodic[1] )
1020        {
1021          out->FaceElements->elementDistribution->numInternal += 2*(numElementsLocal-1-periodic[0]);
1022          out->FaceElements->elementDistribution->numBoundary += 2*(1 + periodic[0]);
1023        }
1024      }
1025      
1026      out->FaceElements->elementDistribution->numLocal = out->FaceElements->elementDistribution->numInternal + out->FaceElements->elementDistribution->numBoundary;
1027    
1028      if( out->FaceElements->elementDistribution->numLocal!=faceNECount )
1029        printf( "ERROR : face element numbers generated %d + %d = %d != %d\n",out->FaceElements->elementDistribution->numInternal, out->FaceElements->elementDistribution->numBoundary, out->FaceElements->elementDistribution->numLocal, faceNECount );
1030    
1031      
1032    
1033      /*  face elements done: */
1034      
1035      /*   condense the nodes: */
1036      Finley_Mesh_resolveNodeIds( out );
1037    
1038      /* prepare mesh for further calculatuions:*/
1039    
1040      //Finley_Mesh_prepare(out) ;
1041    
1042      #ifdef Finley_TRACE
1043      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1044      #endif
1045    
1046      if ( !Finley_MPI_noError( mpi_info )) {
1047          Finley_Mesh_dealloc(out);
1048          return NULL;
1049      }
1050      Paso_MPIInfo_dealloc( mpi_info );
1051    
1052      return out;
1053    }
1054    #endif
1055    
1056    
1057    

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

  ViewVC Help
Powered by ViewVC 1.1.26