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

Diff of /trunk/finley/src/Mesh_rec8.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 26  Line 26 
26    
27  #include "RectangularMesh.h"  #include "RectangularMesh.h"
28    
29    #ifdef PASO_MPI
30    /* get the number of nodes/elements for domain with rank=rank, of size processors
31       where n is the total number of nodes/elements in the global domain */
32    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
33    {
34      rank = size-rank-1;
35    
36      if( rank < n%size )
37        return (index_t)floor(n/size)+1;
38      return (index_t)floor(n/size);
39    }
40    
41    
42    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
43    /* A bit messy, but it only has to be done once... */
44    static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal, dim_t *DOFBoundary )
45    {
46      index_t i0;
47      dim_t numNodesGlobal = numElementsGlobal+1;
48    
49      (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
50      if( rank<size-1 ) // add on node for right hand boundary
51        (*numNodesLocal) += 1;
52    
53      numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
54      periodicLocal[0] = periodicLocal[1] = FALSE;
55      nodesExternal[0] = nodesExternal[1] = 1;
56      if( periodic )
57      {
58        if( size==1 )
59        {
60          numElementsLocal[0] = numElementsGlobal;
61          nodesExternal[0] = nodesExternal[1] = 0;
62          periodicLocal[0] = periodicLocal[1] = TRUE;
63        }
64        else
65        {
66          if( rank==0 )
67          {
68            periodicLocal[0] = TRUE;
69            numNodesLocal[0]++;
70          }
71          if( rank==(size-1) )
72          {      
73            periodicLocal[1] = TRUE;
74            numNodesLocal[0]--;
75            numElementsLocal[0]--;
76          }
77        }
78      }
79      else if( !periodic )
80      {
81        if( rank==0 ){
82          nodesExternal[0]--;
83          numElementsLocal[0]--;
84        }
85        if( rank==(size-1) )
86        {
87          nodesExternal[1]--;
88          numElementsLocal[0]--;
89        }
90      }
91      nodesExternal[0]*=2;
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      firstNode[0] *= 2;
106      
107      numDOFLocal[0] = numNodesLocal[0];
108      if( periodicLocal[0] )
109        numDOFLocal[0]--;
110      
111      DOFExternal[0] = nodesExternal[0];
112      DOFExternal[1] = nodesExternal[1];
113    
114      if( size>1 )
115      {
116        DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
117        DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
118      }
119      else
120      {
121        DOFBoundary[0] = DOFBoundary[1] = 0;
122      }
123    }
124    
125    
126    #endif
127  /**************************************************************/  /**************************************************************/
128    
129  Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {  #ifdef PASO_MPI
130    Finley_Mesh* Finley_RectangularMesh_Rec8_singleCPU(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace, Paso_MPIInfo *mpi_info)
131    #else
132    Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace)
133    #endif
134    {
135    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;
136    index_t k,node0;    index_t k,node0;
137    Finley_Mesh* out=NULL;    Finley_Mesh* out=NULL;
# Line 39  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 142  Finley_Mesh* Finley_RectangularMesh_Rec8
142    N0=2*NE0+1;    N0=2*NE0+1;
143    N1=2*NE1+1;    N1=2*NE1+1;
144    
145    #ifndef PASO_MPI
146    if (N0<=N1) {    if (N0<=N1) {
147       M0=1;       M0=1;
148       M1=N0;       M1=N0;
# Line 46  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 150  Finley_Mesh* Finley_RectangularMesh_Rec8
150       M0=N1;       M0=N1;
151       M1=1;       M1=1;
152    }    }
153    #else
154         M0=1;
155         M1=N0;
156    #endif
157    
158    NFaceElements=0;    NFaceElements=0;
159    if (!periodic[0]) {    if (!periodic[0]) {
# Line 64  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 172  Finley_Mesh* Finley_RectangularMesh_Rec8
172    /*  allocate mesh: */    /*  allocate mesh: */
173        
174    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
175      /* TEMPFIX */  #ifdef PASO_MPI
176  #ifndef PASO_MPI    out=Finley_Mesh_alloc(name,2,order,mpi_info);
177    #else
178    out=Finley_Mesh_alloc(name,2,order);    out=Finley_Mesh_alloc(name,2,order);
179    #endif
180    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
181    
182    #ifdef PASO_MPI
183      out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
184      if (useElementsOnFace) {
185         out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
186         out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
187      } else {
188         out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
189         out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
190      }
191      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
192    #else
193    out->Elements=Finley_ElementFile_alloc(Rec8,out->order);    out->Elements=Finley_ElementFile_alloc(Rec8,out->order);
194    if (useElementsOnFace) {    if (useElementsOnFace) {
195       out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order);       out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order);
# Line 79  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 199  Finley_Mesh* Finley_RectangularMesh_Rec8
199       out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order);       out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order);
200    }    }
201    out->Points=Finley_ElementFile_alloc(Point1,out->order);    out->Points=Finley_ElementFile_alloc(Point1,out->order);
 #else  
   /* TODO */  
202  #endif  #endif
203    
204    if (! Finley_noError()) {    if (! Finley_noError()) {
205         Finley_Mesh_dealloc(out);         Finley_Mesh_dealloc(out);
206         return NULL;         return NULL;
207    }    }
208        
209    /*  allocate tables: */    /*  allocate tables: */
 #ifndef PASO_MPI  
210    Finley_NodeFile_allocTable(out->Nodes,N0*N1);    Finley_NodeFile_allocTable(out->Nodes,N0*N1);
211  #else  #ifdef PASO_MPI
212    /* TODO */    Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1, 0, 0 );
213  #endif  #endif
214    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1);
215    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
# Line 99  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 217  Finley_Mesh* Finley_RectangularMesh_Rec8
217        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
218        return NULL;        return NULL;
219    }    }
220    
221    /*  set nodes: */    /*  set nodes: */
                                                                                                                                                                                                       
222    #pragma omp parallel for private(i0,i1,k)    #pragma omp parallel for private(i0,i1,k)
223    for (i1=0;i1<N1;i1++) {    for (i1=0;i1<N1;i1++) {
224      for (i0=0;i0<N0;i0++) {      for (i0=0;i0<N0;i0++) {
# Line 110  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 228  Finley_Mesh* Finley_RectangularMesh_Rec8
228        out->Nodes->Id[k]=i0+N0*i1;        out->Nodes->Id[k]=i0+N0*i1;
229        out->Nodes->Tag[k]=0;        out->Nodes->Tag[k]=0;
230        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);
231    #ifdef PASO_MPI
232          out->Nodes->Dom[k]=NODE_INTERNAL;
233    #endif
234      }      }
235    }    }
236    /* tags for the faces: */    /* tags for the faces: */
# Line 137  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 258  Finley_Mesh* Finley_RectangularMesh_Rec8
258        out->Elements->Id[k]=k;        out->Elements->Id[k]=k;
259        out->Elements->Tag[k]=0;        out->Elements->Tag[k]=0;
260        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);        out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
261    #ifdef PASO_MPI
262          out->Elements->Dom[k]=ELEMENT_INTERNAL;
263    #endif
264    
265        out->Elements->Nodes[INDEX2(0,k,8)]=node0;        out->Elements->Nodes[INDEX2(0,k,8)]=node0;
266        out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;        out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
# Line 159  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 283  Finley_Mesh* Finley_RectangularMesh_Rec8
283       NUMNODES=3;       NUMNODES=3;
284    }    }
285        
     
286    totalNECount=NE0*NE1;    totalNECount=NE0*NE1;
287    faceNECount=0;    faceNECount=0;
288        
# Line 174  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 297  Finley_Mesh* Finley_RectangularMesh_Rec8
297         out->FaceElements->Id[k]=i0+totalNECount;         out->FaceElements->Id[k]=i0+totalNECount;
298         out->FaceElements->Tag[k]=10;         out->FaceElements->Tag[k]=10;
299         out->FaceElements->Color[k]=i0%2;         out->FaceElements->Color[k]=i0%2;
300    #ifdef PASO_MPI
301             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
302    #endif
303    
304         if (useElementsOnFace) {         if (useElementsOnFace) {
305            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 203  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 329  Finley_Mesh* Finley_RectangularMesh_Rec8
329         out->FaceElements->Id[k]=i0+totalNECount;         out->FaceElements->Id[k]=i0+totalNECount;
330         out->FaceElements->Tag[k]=20;         out->FaceElements->Tag[k]=20;
331         out->FaceElements->Color[k]=i0%2+2;         out->FaceElements->Color[k]=i0%2+2;
332    #ifdef PASO_MPI
333             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
334    #endif
335         if (useElementsOnFace) {         if (useElementsOnFace) {
336            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;
337            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;
# Line 232  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 361  Finley_Mesh* Finley_RectangularMesh_Rec8
361         out->FaceElements->Id[k]=i1+totalNECount;         out->FaceElements->Id[k]=i1+totalNECount;
362         out->FaceElements->Tag[k]=1;         out->FaceElements->Tag[k]=1;
363         out->FaceElements->Color[k]=(i1%2)+4;         out->FaceElements->Color[k]=(i1%2)+4;
364    #ifdef PASO_MPI
365             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
366    #endif
367    
368         if (useElementsOnFace) {         if (useElementsOnFace) {
369            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
# Line 262  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 394  Finley_Mesh* Finley_RectangularMesh_Rec8
394         out->FaceElements->Id[k]=i1+totalNECount;         out->FaceElements->Id[k]=i1+totalNECount;
395         out->FaceElements->Tag[k]=2;         out->FaceElements->Tag[k]=2;
396         out->FaceElements->Color[k]=(i1%2)+6;         out->FaceElements->Color[k]=(i1%2)+6;
397    #ifdef PASO_MPI
398             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
399    #endif
400        
401         if (useElementsOnFace) {         if (useElementsOnFace) {
402            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
# Line 285  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 420  Finley_Mesh* Finley_RectangularMesh_Rec8
420    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
421    out->FaceElements->maxColor=7;    out->FaceElements->maxColor=7;
422    
423    /*  face elements done: */  #ifdef PASO_MPI
424          Finley_ElementFile_setDomainFlags( out->Elements );
425        Finley_ElementFile_setDomainFlags( out->FaceElements );
426        Finley_ElementFile_setDomainFlags( out->ContactElements );
427        Finley_ElementFile_setDomainFlags( out->Points );
428    
429        /* reorder the degrees of freedom */
430        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
431    #endif
432    
433    /*   condense the nodes: */    /*   condense the nodes: */
     
434    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
435        
436    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
   
437    Finley_Mesh_prepare(out) ;    Finley_Mesh_prepare(out) ;
438    
439    
440    #ifdef Finley_TRACE    #ifdef Finley_TRACE
441    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
442    #endif    #endif
# Line 305  Finley_Mesh* Finley_RectangularMesh_Rec8 Line 447  Finley_Mesh* Finley_RectangularMesh_Rec8
447    }    }
448    return out;    return out;
449  }  }
450    
451    #ifdef PASO_MPI
452    Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order,int useElementsOnFace) {
453      dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
454        dim_t N0t, NDOF0t;
455      index_t *numForward=NULL, *numBackward=NULL;
456      index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;
457      Finley_Mesh* out=NULL;
458      char name[50];
459        index_t *indexForward=NULL, *indexBackward=NULL;
460      double time0=Finley_timer();
461      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
462      Paso_MPIInfo *mpi_info = NULL;
463    
464        /* these are used in constructing the face elements, and give the permutation
465           of the node order of a normal element for each face */
466        index_t *facePerm;  
467        index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};
468        index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};
469        index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};
470        index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};
471        index_t face0nc[]={0, 1, 4};
472        index_t face1nc[]={2, 3, 6};
473        index_t face2nc[]={3, 0, 7};
474        index_t face3nc[]={1, 2, 5};
475    
476      NE0=MAX(1,numElements[0]);
477      NE1=MAX(1,numElements[1]);
478      N0=2*NE0+1;
479      N1=2*NE1+1;
480    
481      /* get MPI information */
482      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
483      if (! Finley_noError()) {
484            Finley_Mesh_dealloc(out);
485            return NULL;
486      }
487    
488        // use the serial method for the single CPU case
489        if( mpi_info->size==1 ){
490            out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order,useElementsOnFace, mpi_info);
491            return out;
492        }
493        
494      if( mpi_info->rank==0 )
495        domLeft = TRUE;
496      if( mpi_info->rank==mpi_info->size-1 )
497        domRight = TRUE;
498      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
499        domInternal = TRUE;
500    
501      /* dimensions of the local subdomain */
502      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );  
503    
504      NFaceElements=0;
505      if (!periodic[0]) {
506          NDOF0=N0;
507          NFaceElements+=(domLeft + domRight)*NE1;
508      } else {
509          NDOF0=N0-1;
510      }
511      if (!periodic[1]) {
512          NDOF1=N1;
513          NFaceElements+=2*numElementsLocal;
514      } else {
515          NDOF1=N1-1;
516      }
517    
518        boundaryLeft = !domLeft || periodicLocal[0];
519        boundaryRight = !domRight || periodicLocal[1];
520        N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
521        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
522        firstNodeConstruct = firstNode - 2*boundaryLeft;
523        firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
524    
525      /*  allocate mesh: */
526      
527      sprintf(name,"Rectangular %d x %d mesh",N0,N1);
528      out=Finley_Mesh_alloc(name,2,order,mpi_info);
529    
530      if (! Finley_noError()) return NULL;
531    
532      out->Elements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
533      if (useElementsOnFace) {
534         out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order,mpi_info);
535         out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order,mpi_info);
536      } else {
537         out->FaceElements=Finley_ElementFile_alloc(Line3,out->order,mpi_info);
538         out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order,mpi_info);
539      }
540      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
541    
542      if (! Finley_noError()) {
543           Finley_Mesh_dealloc(out);
544           return NULL;
545      }
546      
547      /*  allocate tables: */
548      Finley_NodeFile_allocTable(out->Nodes,N0t*N1);
549      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );
550      Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);
551      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
552      if (! Finley_noError()) {
553          Finley_Mesh_dealloc(out);
554          return NULL;
555      }
556    
557      /*  set nodes: */
558      #pragma omp parallel for private(i0,i1,k)
559      for (k=0,i1=0;i1<N1;i1++) {
560        for (i0=0;i0<N0t;i0++, k++) {
561          out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];
562          out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];
563          out->Nodes->Id[k]=i0 + i1*N0t;
564          out->Nodes->Tag[k]=0;
565          out->Nodes->Dom[k]=NODE_INTERNAL;
566          out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;
567        }
568            /* take care to make the DOF reflect the internal/external position of the DOF */
569            if( boundaryLeft ){
570                out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;    
571                out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;  
572                out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;  
573                if( periodicLocal[0] ){
574                    out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
575                    out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;  
576                }
577            }
578            if( boundaryRight ){
579                out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;  
580                out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;  
581                out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;  
582            }
583      }
584      /* tags for the faces: */
585      if (!periodic[1]) {
586        for (i0=0;i0<N0t;i0++) {
587          out->Nodes->Tag[i0]+=10;
588          out->Nodes->Tag[i0+N0t*(N1-1)]+=20;
589        }
590      }
591      if (domLeft && !periodic[0]) {
592        for (i1=0;i1<N1;i1++) {
593          out->Nodes->Tag[i1*N0t]+=1;
594        }
595      }
596        if (domRight && !periodic[0]){
597        for (i1=0;i1<N1;i1++) {
598          out->Nodes->Tag[(i1+1)*N0t-1]+=2;
599            }
600        }
601      
602        /* setup the receive information */
603        indexBackward = TMPMEMALLOC( NDOF1*2, index_t );
604        indexForward  = TMPMEMALLOC( NDOF1*2, index_t );
605        if( !(mpi_info->size==2 && periodicLocal[0])){
606            /* Backward (receive) DOF */
607            if( boundaryLeft ){
608                targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;
609    
610                for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
611                    indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
612                    indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
613                    indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
614                }
615    
616                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
617          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
618            }
619            
620            if( boundaryRight ){
621                targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
622    
623                for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
624                    indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
625                    indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
626                    indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
627                }
628    
629                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
630          Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
631            }
632        }
633        else{
634            targetDomain = 1;
635    
636            for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
637                    indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];
638                    indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];
639                    indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];
640            }
641    
642            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );
643            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );
644    
645            for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){
646                    indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];
647                    indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];
648                    indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];
649            }
650    
651            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );
652            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );
653        }
654        TMPMEMFREE( indexBackward );
655        TMPMEMFREE( indexForward );
656    
657      /*   set the elements: */
658      #pragma omp parallel for private(i0,i1,k,node0)
659      for (i1=0;i1<NE1;i1++) {
660        for (i0=0;i0<numElementsLocal;i0++) {
661          k=i0+numElementsLocal*i1;
662                node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t :  2*i0+2*i1*N0t+periodicLocal[0] ;
663    
664          out->Elements->Id[k]=k;
665          out->Elements->Tag[k]=0;
666          out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);
667                out->Elements->Dom[k] = ELEMENT_INTERNAL;
668    
669          out->Elements->Nodes[INDEX2(0,k,8)]=node0;
670          out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;
671          out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;
672          out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;
673          out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;
674          out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;
675          out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;
676          out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;
677        }
678            if( boundaryLeft )
679                out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;
680            if( boundaryRight )
681                out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;
682      }
683      out->Elements->minColor=0;
684      out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);
685      
686        /* set the distribution information for the elements */
687        Finley_ElementFile_setDomainFlags( out->Elements );
688    
689      /*   face elements: */
690      if (useElementsOnFace) {
691         NUMNODES=8;
692      } else {
693         NUMNODES=3;
694      }
695      
696      totalNECount=numElementsLocal*NE1;
697      faceNECount=0;
698        faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);
699      
700      if (!periodic[1]) {
701         /* **  elements on boundary 010 (x2=0): */
702      
703         #pragma omp parallel for private(i0,k,node0)
704         for (i0=0;i0<numElementsLocal; i0++) {
705           k=i0+faceNECount;
706                 kk=i0;
707                 facePerm = useElementsOnFace ? face0 : face0nc;
708    
709           out->FaceElements->Id[k]=i0+totalNECount;
710           out->FaceElements->Tag[k]=10;
711           out->FaceElements->Color[k]=i0%2;
712                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
713                    
714                 for( j=0; j<NUMNODES; j++ )
715                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
716         }
717             if( boundaryLeft ){
718                 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
719                 if( periodicLocal[0] )
720                     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
721             }
722             if( boundaryRight )
723                 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
724         totalNECount+=numElementsLocal;
725         faceNECount+=numElementsLocal;
726      
727         /* **  elements on boundary 020 (x2=1): */
728      
729         #pragma omp parallel for private(i0,k,node0)
730         for (i0=0;i0<numElementsLocal;i0++) {
731           k=i0+faceNECount;
732                 kk=i0 + numElementsLocal*(NE1-1);
733                 facePerm = useElementsOnFace ? face1 : face1nc;
734    
735           out->FaceElements->Id[k]=i0+totalNECount;
736           out->FaceElements->Tag[k]=20;
737           out->FaceElements->Color[k]=i0%2+2;
738                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
739    
740                 for( j=0; j<NUMNODES; j++ )
741                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
742         }
743             if( boundaryLeft ){
744                 out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;
745                 if( periodicLocal[0] )
746                     out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;
747             }
748             if( boundaryRight )
749                 out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;
750         totalNECount+=numElementsLocal;
751         faceNECount+=numElementsLocal;
752      }
753      if (domLeft && !periodicLocal[0]) {
754         /* **  elements on boundary 001 (x1=0): */
755      
756         #pragma omp parallel for private(i1,k,node0)
757         for (i1=0;i1<NE1;i1++) {
758           k=i1+faceNECount;
759                 kk=i1*numElementsLocal;
760                 facePerm = useElementsOnFace ? face2 : face2nc;
761    
762           out->FaceElements->Id[k]=i1+totalNECount;
763           out->FaceElements->Tag[k]=1;
764           out->FaceElements->Color[k]=(i1%2)+4;
765                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
766    
767                 for( j=0; j<NUMNODES; j++ )
768                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
769         }
770         totalNECount+=NE1;
771         faceNECount+=NE1;
772      
773         /* **  elements on boundary 002 (x1=1): */
774      }
775      if (domRight && !periodicLocal[1]) {
776         #pragma omp parallel for private(i1,k,node0)
777         for (i1=0;i1<NE1;i1++) {
778           k=i1+faceNECount;
779                 kk=(i1+1)*numElementsLocal-1;
780                 facePerm = useElementsOnFace ? face3 : face3nc;
781      
782           out->FaceElements->Id[k]=i1+totalNECount;
783           out->FaceElements->Tag[k]=2;
784           out->FaceElements->Color[k]=(i1%2)+6;
785                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
786      
787                 for( j=0; j<NUMNODES; j++ )
788                      out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];
789         }
790         totalNECount+=NE1;
791         faceNECount+=NE1;
792      }
793      out->FaceElements->minColor=0;
794      out->FaceElements->maxColor=7;
795        out->FaceElements->numElements = faceNECount;
796    
797        Finley_ElementFile_setDomainFlags( out->FaceElements );
798    
799      /* setup distribution info for other elements */
800        Finley_ElementFile_setDomainFlags( out->ContactElements );
801        Finley_ElementFile_setDomainFlags( out->Points );
802    
803        /* reorder the degrees of freedom */
804        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
805        
806      /*   condense the nodes: */
807      Finley_Mesh_resolveNodeIds(out);
808      if ( !Finley_MPI_noError( mpi_info )) {
809        Paso_MPIInfo_dealloc( mpi_info );
810        Finley_Mesh_dealloc(out);
811        return NULL;
812      }
813    
814      /* prepare mesh for further calculatuions:*/
815      Finley_Mesh_prepare(out) ;
816    
817      #ifdef Finley_TRACE
818      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
819      #endif
820    
821      if ( !Finley_MPI_noError( mpi_info )) {
822        Paso_MPIInfo_dealloc( mpi_info );
823        Finley_Mesh_dealloc(out);
824        return NULL;
825      }
826    
827        return out;
828    }
829    #endif

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

  ViewVC Help
Powered by ViewVC 1.1.26