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

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

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

revision 730 by bcumming, Mon May 15 04:03:49 2006 UTC revision 782 by bcumming, Tue Jul 18 00:47:47 2006 UTC
# Line 29  Line 29 
29    
30  /**************************************************************/  /**************************************************************/
31    
32  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace) {  #ifdef PASO_MPI
33    /* get the number of nodes/elements for domain with rank=rank, of size processors
34       where n is the total number of nodes/elements in the global domain */
35    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
36    {
37      rank = size-rank-1;
38    
39      if( rank < n%size )
40        return (index_t)floor(n/size)+1;
41      return (index_t)floor(n/size);
42    }
43    
44    
45    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
46    /* A bit messy, but it only has to be done once... */
47    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 )
48    {
49      index_t i0;
50      dim_t numNodesGlobal = numElementsGlobal+1;
51    
52      (*numNodesLocal) = domain_MODdim( rank, size, numNodesGlobal );
53      
54      numElementsLocal[0] = numNodesLocal[0]+1;
55      periodicLocal[0] = periodicLocal[1] = FALSE;
56      nodesExternal[0] = nodesExternal[1] = 1;
57      if( periodic )
58      {
59        if( size==1 )
60        {
61          numElementsLocal[0] = numElementsGlobal;
62          nodesExternal[0] = nodesExternal[1] = 0;
63          periodicLocal[0] = periodicLocal[1] = TRUE;
64        }
65        else
66        {
67          if( rank==0 )
68          {
69            periodicLocal[0] = TRUE;
70            numNodesLocal[0]++;
71          }
72          if( rank==(size-1) )
73          {      
74            periodicLocal[1] = TRUE;
75            numNodesLocal[0]--;
76            numElementsLocal[0]--;
77          }
78        }
79      }
80      else if( !periodic )
81      {
82        if( rank==0 ){
83          nodesExternal[0]--;
84          numElementsLocal[0]--;
85        }
86        if( rank==(size-1) )
87        {
88          nodesExternal[1]--;
89          numElementsLocal[0]--;
90        }
91      }
92      numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
93      
94      numElementsInternal[0] = numElementsLocal[0];
95      if( (rank==0) && (rank==size-1) );
96        
97      else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
98          numElementsInternal[0] -= 1;
99      else
100        numElementsInternal[0] -= 2;
101    
102      firstNode[0] = 0;
103      for( i0=0; i0<rank; i0++ )
104        firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
105    
106      numDOFLocal[0] = numNodesLocal[0];
107      if( periodicLocal[0] )
108      {
109        numDOFLocal[0]--;
110      }
111      DOFExternal[0] = nodesExternal[0];
112      DOFExternal[1] = nodesExternal[1];
113    }
114    
115    #endif
116    
117    #ifdef PASO_MPI
118    Finley_Mesh* Finley_RectangularMesh_Hex8_singleCPU(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace,Paso_MPIInfo *mpi_info)
119    #else
120    Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
121    #endif
122    {
123    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
124    index_t node0;    index_t node0;
125    Finley_Mesh* out;    Finley_Mesh* out;
# Line 98  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 188  Finley_Mesh* Finley_RectangularMesh_Hex8
188    /*  allocate mesh: */    /*  allocate mesh: */
189        
190    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
191    /* TEMPFIX */  
192  #ifndef PASO_MPI  #ifndef PASO_MPI
193    out=Finley_Mesh_alloc(name,3,order);    out=Finley_Mesh_alloc(name,3,order);
194  #else  #else
195    out=Finley_Mesh_alloc(name,3,order,NULL);    out=Finley_Mesh_alloc(name,3,order,mpi_info);
196  #endif  #endif
197    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
198    
199  #ifndef PASO_MPI  #ifdef PASO_MPI
200      out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
201      if (useElementsOnFace) {
202         out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
203         out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
204      } else {
205         out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
206         out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
207      }
208      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
209    #else
210    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);
211    if (useElementsOnFace) {    if (useElementsOnFace) {
212       out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);       out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);
# Line 116  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 216  Finley_Mesh* Finley_RectangularMesh_Hex8
216       out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);       out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);
217    }    }
218    out->Points=Finley_ElementFile_alloc(Point1,out->order);    out->Points=Finley_ElementFile_alloc(Point1,out->order);
 #else  
   out->Elements=Finley_ElementFile_alloc(Hex8,out->order,NULL);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,NULL);  
      out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,NULL);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,NULL);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,NULL);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order,NULL);  
219  #endif  #endif
220    if (! Finley_noError()) {    if (! Finley_noError()) {
221        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
# Line 133  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 223  Finley_Mesh* Finley_RectangularMesh_Hex8
223    }    }
224    
225        
226    /*  allocate tables: */    /*  allocate tables: */
 #ifndef PASO_MPI    
227    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
228  #else  #ifdef PASO_MPI
229   /* TODO */    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
230  #endif  #endif
231    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
232    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
# Line 159  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 248  Finley_Mesh* Finley_RectangularMesh_Hex8
248          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
249          out->Nodes->Tag[k]=0;          out->Nodes->Tag[k]=0;
250          out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);          out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
251    #ifdef PASO_MPI
252            out->Nodes->Dom[k]=NODE_INTERNAL;
253    #endif
254        }        }
255      }      }
256    }    }
# Line 201  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 293  Finley_Mesh* Finley_RectangularMesh_Hex8
293          out->Elements->Id[k]=k;          out->Elements->Id[k]=k;
294          out->Elements->Tag[k]=0;          out->Elements->Tag[k]=0;
295          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
296    #ifdef PASO_MPI
297            out->Elements->Dom[k]=ELEMENT_INTERNAL;
298    #endif
299    
300          out->Elements->Nodes[INDEX2(0,k,8)]=node0;          out->Elements->Nodes[INDEX2(0,k,8)]=node0;
301          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;          out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
# Line 240  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 335  Finley_Mesh* Finley_RectangularMesh_Hex8
335           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
336           out->FaceElements->Tag[k]=100;           out->FaceElements->Tag[k]=100;
337           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);           out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
338    #ifdef PASO_MPI
339            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
340    #endif
341    
342           if  (useElementsOnFace) {           if  (useElementsOnFace) {
343              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 272  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 370  Finley_Mesh* Finley_RectangularMesh_Hex8
370           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
371           out->FaceElements->Tag[k]=200;           out->FaceElements->Tag[k]=200;
372           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
373    #ifdef PASO_MPI
374            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
375    #endif
376    
377           if  (useElementsOnFace) {           if  (useElementsOnFace) {
378              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;
# Line 307  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 408  Finley_Mesh* Finley_RectangularMesh_Hex8
408           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
409           out->FaceElements->Tag[k]=1;           out->FaceElements->Tag[k]=1;
410           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
411    #ifdef PASO_MPI
412            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
413    #endif
414    
415           if  (useElementsOnFace) {           if  (useElementsOnFace) {
416              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 339  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 443  Finley_Mesh* Finley_RectangularMesh_Hex8
443           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
444           out->FaceElements->Tag[k]=2;           out->FaceElements->Tag[k]=2;
445           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
446    #ifdef PASO_MPI
447            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
448    #endif
449    
450           if  (useElementsOnFace) {           if  (useElementsOnFace) {
451              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;
# Line 372  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 479  Finley_Mesh* Finley_RectangularMesh_Hex8
479           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
480           out->FaceElements->Tag[k]=10;           out->FaceElements->Tag[k]=10;
481           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
482    #ifdef PASO_MPI
483            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
484    #endif
485    
486           if  (useElementsOnFace) {           if  (useElementsOnFace) {
487              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 404  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 514  Finley_Mesh* Finley_RectangularMesh_Hex8
514           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
515           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
516           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
517    #ifdef PASO_MPI
518            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
519    #endif
520    
521           if  (useElementsOnFace) {           if  (useElementsOnFace) {
522              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;
# Line 429  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 542  Finley_Mesh* Finley_RectangularMesh_Hex8
542    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
543    out->FaceElements->maxColor=23;    out->FaceElements->maxColor=23;
544        
545    /*  face elements done: */  #ifdef PASO_MPI
546          Finley_ElementFile_setDomainFlags( out->Elements );
547        Finley_ElementFile_setDomainFlags( out->FaceElements );
548        Finley_ElementFile_setDomainFlags( out->ContactElements );
549        Finley_ElementFile_setDomainFlags( out->Points );
550    
551        /* reorder the degrees of freedom */
552        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
553    #endif
554        
555    /*   condense the nodes: */    /*   condense the nodes: */
     
556    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
557    
558    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
# Line 448  Finley_Mesh* Finley_RectangularMesh_Hex8 Line 568  Finley_Mesh* Finley_RectangularMesh_Hex8
568    }    }
569    return out;    return out;
570  }  }
571    
572    #ifdef PASO_MPI
573    Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,double* Length,bool_t* periodic, index_t order,bool_t useElementsOnFace)
574    {
575      dim_t N0,N1,N2,N0t,NDOF0t,NE0,NE1,NE2,i0,i1,i2,kk,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;//,M0,M1,M2;
576      dim_t idCount, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal;
577      bool_t dom_left, dom_right, dom_internal;
578      index_t firstNode=0, DOFcount=0, node0, node1, node2;
579      index_t targetDomain=-1, firstNodeConstruct, j;
580      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
581        index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
582      Finley_Mesh* out;
583    
584      char name[50];
585      Paso_MPIInfo *mpi_info = NULL;
586      double time0=Finley_timer();
587    
588        index_t face0[] = {3, 0, 4, 7, 2, 1, 5, 6};
589        index_t face1[] = {1, 2, 6, 5, 0, 3, 7, 4};
590        index_t face2[] = {0, 3, 2, 1, 4, 7, 6, 5};
591        index_t face3[] = {4, 5, 6, 7, 0, 1, 2, 3};
592        index_t face4[] = {0, 1, 5, 4, 3, 2, 6, 7};
593        index_t face5[] = {3, 7, 6, 2, 0, 4, 5, 1};
594      NE0=MAX(1,numElements[0]);
595      NE1=MAX(1,numElements[1]);
596      NE2=MAX(1,numElements[2]);
597      N0=NE0+1;
598      N1=NE1+1;
599      N2=NE2+1;
600    
601    
602      /* get MPI information */
603      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
604      if (! Finley_noError())
605            return NULL;
606    
607      /* use the serial version to generate the mesh for the 1-CPU case */
608      if( mpi_info->size==1 )
609      {
610        out =  Finley_RectangularMesh_Hex8_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info );
611            return out;
612      }    
613    
614      if( mpi_info->rank==0 )
615        domLeft = TRUE;
616      if( mpi_info->rank==mpi_info->size-1 )
617        domRight = TRUE;
618      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
619        domInternal = TRUE;
620    
621      /* dimensions of the local subdomain */
622      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal );  
623    
624      /* count Degrees of Freedom along each axis */
625      NDOF0 = N0 - periodic[0];
626      NDOF1 = N1 - periodic[1];
627      NDOF2 = N2 - periodic[2];
628    
629      /* count face elements */
630      /* internal face elements */
631      NFaceElements = 0;
632      if( !periodic[0] )
633        NFaceElements += (domLeft+domRight)*NE1*NE2;
634      if( !periodic[1] )
635        NFaceElements += 2*numElementsLocal*NE2;
636      if( !periodic[2] )
637        NFaceElements += 2*numElementsLocal*NE1;
638      
639        boundaryLeft = !domLeft || periodicLocal[0];
640        boundaryRight = !domRight || periodicLocal[1];
641        N0t = numNodesLocal + boundaryRight + boundaryLeft;
642        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft;
643        firstNodeConstruct = firstNode - boundaryLeft;
644        firstNodeConstruct = firstNodeConstruct<0 ? N0-2 : firstNodeConstruct;
645    
646      /*  allocate mesh: */
647      sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
648    
649      out=Finley_Mesh_alloc(name,3,order,mpi_info);
650      if (! Finley_noError()) return NULL;
651    
652      out->Elements=Finley_ElementFile_alloc(Hex8,out->order,mpi_info);
653      if (useElementsOnFace) {
654         out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order,mpi_info);
655         out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order,mpi_info);
656      } else {
657         out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order,mpi_info);
658         out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order,mpi_info);
659      }
660      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
661    
662      if (! Finley_noError()) {
663          Finley_Mesh_dealloc(out);
664          return NULL;
665      }
666      
667      /*  allocate tables: */
668      Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
669      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*2, 0 );
670      Finley_ElementFile_allocTable(out->Elements,(numElementsLocal)*NE1*NE2);
671      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
672      if (! Finley_noError()) {
673          Finley_Mesh_dealloc(out);
674          return NULL;
675      }
676      
677      /*  set nodes: */
678      /* INTERNAL/BOUNDARY NODES */
679        k=0;
680      #pragma omp parallel for private(i0,i1,i2,k)
681      for (i2=0;i2<N2;i2++) {
682        for (i1=0;i1<N1;i1++) {
683          for (i0=0;i0<N0t;i0++,k++) {        
684            out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
685            out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
686            out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
687            out->Nodes->Id[k]=k;
688            out->Nodes->Tag[k]=0;
689            out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
690                    out->Nodes->Dom[k]=NODE_INTERNAL;
691          }
692        }
693      }
694    
695        /* mark the nodes that reference external and boundary DOF as such */
696        if( boundaryLeft ){
697            for( i1=0; i1<N1; i1++ )
698                for( i2=0; i2<N2; i2++ ) {
699                    out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
700                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_BOUNDARY;
701                }
702        }
703        if( boundaryRight ){
704            for( i1=0; i1<N1; i1++ )
705                for( i2=0; i2<N2; i2++ ) {
706                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
707                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
708                }
709        }
710        if( periodicLocal[0] ){
711            for( i1=0; i1<N1; i1++ )
712                for( i2=0; i2<N2; i2++ ) {
713                    out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+2];
714                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
715                }
716        }
717            
718      /* tag Nodes that are referenced by face elements */
719      if (!periodic[2]) {    
720        for (i1=0;i1<N1;i1++) {
721          for (i0=0;i0<N0t;i0++) {  
722             out->Nodes->Tag[i0 + N0t*i1]+=100;
723             out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
724           }
725         }
726      }
727      if (!periodic[1]) {
728        for (i2=0;i2<N2;i2++) {
729          for (i0=0;i0<N0t;i0++) {
730             out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
731             out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
732          }
733        }
734      }
735      if (!periodic[0] && !domInternal ) {
736        for (i2=0;i2<N2;i2++) {
737          for (i1=0;i1<N1;i1++) {
738            if( domLeft )
739              out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
740            if( domRight )
741              out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
742          }
743        }
744      }
745    
746        /* form the boudary communication information */
747        forwardDOF  = MEMALLOC(NDOF1*NDOF2,index_t);
748        backwardDOF = MEMALLOC(NDOF1*NDOF2,index_t);
749        if( !(mpi_info->size==2 && periodicLocal[0])){
750            if( boundaryLeft  ) {
751                targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
752                for( i2=0; i2<NDOF2; i2++ ){
753                    for( i1=0; i1<NDOF1; i1++ ){
754                        forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
755                        backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
756                    }
757                }
758                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
759                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
760            }
761            if( boundaryRight ) {
762                targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
763                for( i2=0; i2<NDOF2; i2++ ){
764                    for( i1=0; i1<NDOF1; i1++ ){
765                        forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
766                        backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
767                    }
768                }
769                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
770                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
771            }
772        } else{
773            /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
774            targetDomain = 1;
775            
776            for( i2=0; i2<NDOF2; i2++ ){
777                for( i1=0; i1<NDOF1; i1++ ){
778                    forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
779                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
780                }
781            }
782            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
783            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
784    
785            for( i2=0; i2<NDOF2; i2++ ){
786                for( i1=0; i1<NDOF1; i1++ ){
787                    forwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
788                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
789                }
790            }
791            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
792            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
793        }
794        MEMFREE( forwardDOF );
795        MEMFREE( backwardDOF );
796      /*   set the elements: */
797    
798      /* INTERNAL elements */
799      k = 0;
800      #pragma omp parallel for private(i0,i1,i2,k,node0)
801      for (i2=0;i2<NE2;i2++) {
802        for (i1=0;i1<NE1;i1++) {
803          for (i0=0;i0<numElementsLocal;i0++,k++) {
804                    node0 = (periodicLocal[0] && !i0) ? i1*N0t + i2*N1*N0t :  i1*N0t + i2*N1*N0t + i0 + periodicLocal[0];
805            out->Elements->Id[k]=k;
806                    
807            out->Elements->Tag[k]=0;
808            out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;
809                    out->Elements->Dom[k]=ELEMENT_INTERNAL;
810    
811            out->Elements->Nodes[INDEX2(0,k,8)]=node0;
812            out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;
813            out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0t+1;
814            out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0t;
815            out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0t*N1;
816            out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t*N1+1;
817            out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0t*N1+N0t+1;
818            out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t*N1+N0t;
819    
820          }
821        }
822      }
823        out->Elements->minColor=0;
824      out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
825        if( boundaryLeft )
826            for( i2=0; i2<NE2; i2++ )
827                for( i1=0; i1<NE1; i1++ )
828                    out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
829        if( boundaryRight )
830            for( i2=0; i2<NE2; i2++ )
831                for( i1=0; i1<NE1; i1++ )
832                    out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
833            
834        Finley_ElementFile_setDomainFlags( out->Elements );
835        
836      /*   face elements: */
837      if  (useElementsOnFace) {
838         NUMNODES=8;
839      } else {
840         NUMNODES=4;
841      }
842      totalNECount=out->Elements->numElements;
843      faceNECount=0;
844        idCount = totalNECount;
845      
846      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
847        numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
848      if (!periodic[2]) {
849         /*  elements on boundary 100 (x3=0): */
850      
851         #pragma omp parallel for private(i0,i1,k)
852         for (i1=0;i1<NE1;i1++) {
853           for (i0=0; i0<numElementsLocal; i0++) {
854             k=i0+numElementsLocal*i1+faceNECount;
855                     kk=i0 + i1*numElementsLocal;
856                     facePerm = face2;
857      
858             out->FaceElements->Id[k]=idCount++;
859             out->FaceElements->Tag[k]=100;
860             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
861             out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
862    
863                     for( j=0; j<NUMNODES; j++ )
864                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
865           }
866         }
867           if( boundaryLeft ){
868                for( i1=0; i1<NE1; i1++ )
869                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
870                if( periodicLocal[0] )
871                    for( i1=0; i1<NE1; i1++ )
872                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
873             }
874           if( boundaryRight )
875                for( i1=0; i1<NE1; i1++ )
876                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
877         totalNECount+=NE1*numElementsLocal;
878         faceNECount+=NE1*numElementsLocal;
879            
880         /* **  elements on boundary 200 (x3=1) */
881      
882         #pragma omp parallel for private(i0,i1,k)
883         for (i1=0;i1<NE1;i1++) {
884           for (i0=0;i0<numElementsLocal;i0++) {
885             k=i0+numElementsLocal*i1+faceNECount;
886                     kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
887                     facePerm = face3;
888      
889             out->FaceElements->Id[k]=idCount++;
890             out->FaceElements->Tag[k]=200;
891             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
892             out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
893    
894                     for( j=0; j<NUMNODES; j++ )
895                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
896           }
897         }
898           if( boundaryLeft ){
899                for( i1=0; i1<NE1; i1++ )
900                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
901                if( periodicLocal[0] )
902                    for( i1=0; i1<NE1; i1++ )
903                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
904             }
905           if( boundaryRight )
906                for( i1=0; i1<NE1; i1++ )
907                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
908         totalNECount+=NE1*numElementsLocal;
909         faceNECount+=NE1*numElementsLocal;
910      }
911      if (!periodic[0] && !domInternal) {
912         /* **  elements on boundary 001 (x1=0): */
913         if( domLeft ){
914                 #pragma omp parallel for private(i1,i2,k)
915                 for (i2=0;i2<NE2;i2++) {
916                     for (i1=0;i1<NE1;i1++) {
917                         k=i1+NE1*i2+faceNECount;
918                         kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
919                         facePerm = face0;
920            
921                         out->FaceElements->Id[k]=idCount++;
922                         out->FaceElements->Tag[k]=1;
923                         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
924                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
925    
926                         for( j=0; j<NUMNODES; j++ )
927                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
928                     }
929                 }
930                 totalNECount+=NE1*NE2;
931                 faceNECount+=NE1*NE2;
932         }
933         /* **  elements on boundary 002 (x1=1): */
934             if( domRight ) {
935                 #pragma omp parallel for private(i1,i2,k)
936                 for (i2=0;i2<NE2;i2++) {
937                     for (i1=0;i1<NE1;i1++) {
938                         k=i1+NE1*i2+faceNECount;
939                         kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
940                         facePerm = face1;
941            
942                         out->FaceElements->Id[k]=idCount++;
943                         out->FaceElements->Tag[k]=2;
944                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
945                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
946    
947                         for( j=0; j<NUMNODES; j++ )
948                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
949                     }
950                 }
951           totalNECount+=NE1*NE2;
952           faceNECount+=NE1*NE2;
953         }
954      }
955      if (!periodic[1]) {
956         /* **  elements on boundary 010 (x2=0): */
957      
958         #pragma omp parallel for private(i0,i2,k)
959         for (i2=0;i2<NE2;i2++) {
960           for (i0=0;i0<numElementsLocal;i0++) {
961             k=i0+numElementsLocal*i2+faceNECount;
962                     kk=i0+numElementsLocal*NE1*i2;
963                     facePerm = face4;
964      
965             out->FaceElements->Id[k]=idCount++;
966             out->FaceElements->Tag[k]=10;
967             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
968             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
969    
970                     for( j=0; j<NUMNODES; j++ )
971                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
972           }
973         }
974           if( boundaryLeft ){
975                for( i2=0; i2<NE2; i2++ )
976                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
977                if( periodicLocal[0] )
978                    for( i2=0; i2<NE2; i2++ )
979                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
980             }
981           if( boundaryRight )
982                for( i2=0; i2<NE2; i2++ )
983                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
984         totalNECount+=numElementsLocal*NE2;
985         faceNECount+=numElementsLocal*NE2;
986      
987         /* **  elements on boundary 020 (x2=1): */
988      
989         #pragma omp parallel for private(i0,i2,k)
990         for (i2=0;i2<NE2;i2++) {
991           for (i0=0;i0<numElementsLocal;i0++) {
992             k=i0+numElementsLocal*i2+faceNECount;
993                     kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
994                     facePerm = face5;
995      
996             out->FaceElements->Tag[k]=20;
997             out->FaceElements->Id[k]=idCount++;
998             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
999             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1000    
1001                     for( j=0; j<NUMNODES; j++ )
1002                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
1003           }
1004         }
1005           if( boundaryLeft ){
1006                for( i2=0; i2<NE2; i2++ )
1007                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1008                if( periodicLocal[0] )
1009                    for( i2=0; i2<NE2; i2++ )
1010                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1011             }
1012           if( boundaryRight )
1013                for( i2=0; i2<NE2; i2++ )
1014                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1015         totalNECount+=numElementsLocal*NE2;
1016         faceNECount+=numElementsLocal*NE2;
1017      }
1018        out->FaceElements->elementDistribution->numInternal = faceNECount;
1019        
1020      out->FaceElements->minColor=0;
1021      out->FaceElements->maxColor=23;
1022        out->FaceElements->numElements=faceNECount;
1023        
1024        Finley_ElementFile_setDomainFlags( out->FaceElements );
1025    
1026      /* setup distribution info for other elements */
1027        Finley_ElementFile_setDomainFlags( out->ContactElements );
1028        Finley_ElementFile_setDomainFlags( out->Points );
1029    
1030        /* reorder the degrees of freedom */
1031        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1032    
1033      /*   condense the nodes: */
1034      Finley_Mesh_resolveNodeIds(out);
1035      if( !Finley_MPI_noError(mpi_info) )
1036      {
1037        Paso_MPIInfo_dealloc( mpi_info );
1038        Finley_Mesh_dealloc(out);
1039        return NULL;
1040      }
1041    
1042      /* prepare mesh for further calculatuions:*/
1043      Finley_Mesh_prepare(out);
1044      if( !Finley_MPI_noError(mpi_info) )
1045      {
1046        Paso_MPIInfo_dealloc( mpi_info );
1047        Finley_Mesh_dealloc(out);
1048        return NULL;
1049      }
1050    
1051      /* free up memory */
1052      Paso_MPIInfo_dealloc( mpi_info );
1053    
1054      #ifdef Finley_TRACE
1055      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1056      #endif
1057    
1058      return out;
1059    }
1060    #endif
1061    

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

  ViewVC Help
Powered by ViewVC 1.1.26