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

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

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

trunk/finley/src/finleyC/Mesh_hex20.c revision 155 by jgs, Wed Nov 9 02:02:19 2005 UTC trunk/finley/src/Mesh_hex20.c revision 787 by bcumming, Wed Jul 26 01:46:45 2006 UTC
# Line 1  Line 1 
1  /*  /*
2   ******************************************************************************   ************************************************************
3   *                                                                            *   *          Copyright 2006 by ACcESS MNRF                   *
4   *       COPYRIGHT  ACcESS 2003,2004,2005 -  All Rights Reserved              *   *                                                          *
5   *                                                                            *   *              http://www.access.edu.au                    *
6   * This software is the property of ACcESS. No part of this code              *   *       Primary Business: Queensland, Australia            *
7   * may be copied in any form or by any means without the expressed written    *   *  Licensed under the Open Software License version 3.0    *
8   * consent of ACcESS.  Copying, use or modification of this software          *   *     http://www.opensource.org/licenses/osl-3.0.php       *
9   * by any unauthorised person is illegal unless that person has a software    *   *                                                          *
10   * license agreement with ACcESS.                                             *   ************************************************************
  *                                                                            *  
  ******************************************************************************  
11  */  */
12    
13  /**************************************************************/  /**************************************************************/
# Line 30  Line 28 
28    
29  #include "RectangularMesh.h"  #include "RectangularMesh.h"
30    
31  /**************************************************************/  #ifdef PASO_MPI
32    /* get the number of nodes/elements for domain with rank=rank, of size processors
33       where n is the total number of nodes/elements in the global domain */
34    static index_t domain_MODdim( index_t rank, index_t size, index_t n )
35    {
36      rank = size-rank-1;
37    
38      if( rank < n%size )
39        return (index_t)floor(n/size)+1;
40      return (index_t)floor(n/size);
41    }
42    
43    
44    /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */
45    /* A bit messy, but it only has to be done once... */
46    static void domain_calculateDimension( index_t rank, dim_t size, dim_t numElementsGlobal, bool_t periodic, dim_t *numNodesLocal, dim_t *numDOFLocal, dim_t *numElementsLocal, dim_t *numElementsInternal, dim_t *firstNode, dim_t *nodesExternal, dim_t *DOFExternal, dim_t *numNodesExternal, bool_t *periodicLocal, dim_t *DOFBoundary )
47    {
48      index_t i0;
49      dim_t numNodesGlobal = numElementsGlobal+1;
50    
51      (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;
52      if( rank<size-1 ) // add on node for right hand boundary
53        (*numNodesLocal) += 1;
54    
55      numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;
56      periodicLocal[0] = periodicLocal[1] = FALSE;
57      nodesExternal[0] = nodesExternal[1] = 1;
58      if( periodic )
59      {
60        if( size==1 )
61        {
62          numElementsLocal[0] = numElementsGlobal;
63          nodesExternal[0] = nodesExternal[1] = 0;
64          periodicLocal[0] = periodicLocal[1] = TRUE;
65        }
66        else
67        {
68          if( rank==0 )
69          {
70            periodicLocal[0] = TRUE;
71            numNodesLocal[0]++;
72          }
73          if( rank==(size-1) )
74          {      
75            periodicLocal[1] = TRUE;
76            numNodesLocal[0]--;
77            numElementsLocal[0]--;
78          }
79        }
80      }
81      else if( !periodic )
82      {
83        if( rank==0 ){
84          nodesExternal[0]--;
85          numElementsLocal[0]--;
86        }
87        if( rank==(size-1) )
88        {
89          nodesExternal[1]--;
90          numElementsLocal[0]--;
91        }
92      }
93      nodesExternal[0]*=2;
94      numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];
95      
96      numElementsInternal[0] = numElementsLocal[0];
97      if( (rank==0) && (rank==size-1) );
98        
99      else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )
100          numElementsInternal[0] -= 1;
101      else
102        numElementsInternal[0] -= 2;
103    
104      firstNode[0] = 0;
105      for( i0=0; i0<rank; i0++ )
106        firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );
107      firstNode[0] *= 2;
108      
109      numDOFLocal[0] = numNodesLocal[0];
110      if( periodicLocal[0] )
111        numDOFLocal[0]--;
112      
113      DOFExternal[0] = nodesExternal[0];
114      DOFExternal[1] = nodesExternal[1];
115    
116      if( size>1 )
117      {
118        DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;
119        DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;
120      }
121      else
122      {
123        DOFBoundary[0] = DOFBoundary[1] = 0;
124      }
125    }
126    #endif
127    /**************************************************************/
128    #ifdef PASO_MPI
129    Finley_Mesh* Finley_RectangularMesh_Hex20_singleCPU(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace, Paso_MPIInfo *mpi_info) {
130    #else
131  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {
132    #endif
133    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;
134    index_t node0;    index_t node0;
135    Finley_Mesh* out;    Finley_Mesh* out;
# Line 101  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 198  Finley_Mesh* Finley_RectangularMesh_Hex2
198    /*  allocate mesh: */    /*  allocate mesh: */
199        
200    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
201    #ifdef PASO_MPI
202      out=Finley_Mesh_alloc(name,3,order,mpi_info);
203    
204      if (! Finley_noError()) return NULL;
205    
206      out->Elements=Finley_ElementFile_alloc(Hex20,out->order,mpi_info);
207      if (useElementsOnFace) {
208         out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,mpi_info);
209         out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,mpi_info);
210      } else {
211         out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
212         out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,mpi_info);
213      }
214      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
215      if (! Finley_noError()) {
216          Finley_Mesh_dealloc(out);
217          return NULL;
218      }
219    
220      
221      /*  allocate tables: */
222      Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
223      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
224      Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
225      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
226    #else
227    out=Finley_Mesh_alloc(name,3,order);    out=Finley_Mesh_alloc(name,3,order);
228    
229    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
230    
231    out->Elements=Finley_ElementFile_alloc(Hex20,out->order);    out->Elements=Finley_ElementFile_alloc(Hex20,out->order);
# Line 120  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 244  Finley_Mesh* Finley_RectangularMesh_Hex2
244    
245        
246    /*  allocate tables: */    /*  allocate tables: */
     
247    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
248    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
249    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
250    #endif
251    if (! Finley_noError()) {    if (! Finley_noError()) {
252        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
253        return NULL;        return NULL;
254    }    }
255    
256        /* create nodes */
257    #pragma omp parallel for private(i0,i1,i2,k)    #pragma omp parallel for private(i0,i1,i2,k)
258    for (i2=0;i2<N2;i2++) {    for (i2=0;i2<N2;i2++) {
259      for (i1=0;i1<N1;i1++) {      for (i1=0;i1<N1;i1++) {
# Line 140  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 265  Finley_Mesh* Finley_RectangularMesh_Hex2
265          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
266          out->Nodes->Tag[k]=0;          out->Nodes->Tag[k]=0;
267          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);
268    #ifdef PASO_MPI
269            out->Nodes->Dom[k]=NODE_INTERNAL;
270    #endif
271        }        }
272      }      }
273    }    }
# Line 183  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 311  Finley_Mesh* Finley_RectangularMesh_Hex2
311          out->Elements->Id[k]=k;          out->Elements->Id[k]=k;
312          out->Elements->Tag[k]=0;          out->Elements->Tag[k]=0;
313          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);
314    #ifdef PASO_MPI
315            out->Elements->Dom[k]=ELEMENT_INTERNAL;
316    #endif
317    
318          out->Elements->Nodes[INDEX2(0,k,20)]=node0;          out->Elements->Nodes[INDEX2(0,k,20)]=node0;
319          out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;          out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
# Line 233  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 364  Finley_Mesh* Finley_RectangularMesh_Hex2
364          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
365          out->FaceElements->Tag[k]=100;          out->FaceElements->Tag[k]=100;
366          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
367    #ifdef PASO_MPI
368            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
369    #endif
370        
371          if  (useElementsOnFace) {          if  (useElementsOnFace) {
372             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 280  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 414  Finley_Mesh* Finley_RectangularMesh_Hex2
414          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
415          out->FaceElements->Tag[k]=200;          out->FaceElements->Tag[k]=200;
416          out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;          out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
417    #ifdef PASO_MPI
418            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
419    #endif
420    
421          if  (useElementsOnFace) {          if  (useElementsOnFace) {
422             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
# Line 334  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 471  Finley_Mesh* Finley_RectangularMesh_Hex2
471           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
472           out->FaceElements->Tag[k]=1;           out->FaceElements->Tag[k]=1;
473           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
474    #ifdef PASO_MPI
475            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
476    #endif
477    
478           if  (useElementsOnFace) {           if  (useElementsOnFace) {
479              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 387  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 527  Finley_Mesh* Finley_RectangularMesh_Hex2
527           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
528           out->FaceElements->Tag[k]=2;           out->FaceElements->Tag[k]=2;
529           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
530    #ifdef PASO_MPI
531            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
532    #endif
533    
534           if  (useElementsOnFace) {           if  (useElementsOnFace) {
535              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
# Line 442  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 585  Finley_Mesh* Finley_RectangularMesh_Hex2
585           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
586           out->FaceElements->Tag[k]=10;           out->FaceElements->Tag[k]=10;
587           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
588    #ifdef PASO_MPI
589            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
590    #endif
591    
592           if  (useElementsOnFace) {           if  (useElementsOnFace) {
593              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 495  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 641  Finley_Mesh* Finley_RectangularMesh_Hex2
641           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
642           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
643           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
644    #ifdef PASO_MPI
645            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
646    #endif
647    
648           if  (useElementsOnFace) {           if  (useElementsOnFace) {
649              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
# Line 539  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 688  Finley_Mesh* Finley_RectangularMesh_Hex2
688    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
689    out->FaceElements->maxColor=24;    out->FaceElements->maxColor=24;
690    
691    /*  face elements done: */  #ifdef PASO_MPI
692          Finley_ElementFile_setDomainFlags( out->Elements );
693        Finley_ElementFile_setDomainFlags( out->FaceElements );
694        Finley_ElementFile_setDomainFlags( out->ContactElements );
695        Finley_ElementFile_setDomainFlags( out->Points );
696    
697        /* reorder the degrees of freedom */
698        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
699    #endif
700        
701    /*   condense the nodes: */    /*   condense the nodes: */
     
702    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
703    
704    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
# Line 558  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 714  Finley_Mesh* Finley_RectangularMesh_Hex2
714    }    }
715    return out;    return out;
716  }  }
717    #ifdef PASO_MPI
718    Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {
719      dim_t N0,N0t,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF0t,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
720      dim_t kk,iI, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
721        bool_t dom_left, dom_right, dom_internal;
722      index_t firstNode=0, DOFcount=0, node0, node1, node2, idCount;
723      index_t targetDomain=-1, firstNodeConstruct, j;
724      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
725        index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
726      Finley_Mesh* out;
727      char name[50];
728      double time0=Finley_timer();
729      Paso_MPIInfo *mpi_info = NULL;
730    
731      NE0=MAX(1,numElements[0]);
732      NE1=MAX(1,numElements[1]);
733      NE2=MAX(1,numElements[2]);
734      N0=2*NE0+1;
735      N1=2*NE1+1;
736      N2=2*NE2+1;
737    
738        index_t face0[] = {0,4,7,3,1,5,6,2,12,19,15,11,8,16,18,10,13,17,14,9};
739        index_t face1[] = {1,2,6,5,0,3,7,4,9,14,17,13,8,10,18,16,11,15,19,12};
740        index_t face2[] = {0,1,5,4,3,2,6,7,8,13,16,12,11,9,17,19,14,18,15,10};
741        index_t face3[] = {3,7,6,2,0,4,5,1,15,18,14,10,11,19,17,9,12,16,13,8};
742        index_t face4[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,19,18,17,16};
743        index_t face5[] = {4,5,6,7,0,1,2,3,16,17,18,19,12,13,14,15,8,9,10,11};
744        
745        index_t face0R[] = {0,4,7,3,12,19,15,11};
746        index_t face1R[] = {1,2,6,5,9,14,17,13};
747        index_t face2R[] = {0,1,5,4,8,13,16,12};
748        index_t face3R[] = {3,7,6,2,15,18,14,10};
749        index_t face4R[] = {0,3,2,1,11,10,9,8};
750        index_t face5R[] = {4,5,6,7,16,17,18,19};
751        
752      /* get MPI information */
753      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
754      if (! Finley_noError())
755            return NULL;
756    
757      /* use the serial version to generate the mesh for the 1-CPU case */
758      if( mpi_info->size==1 )
759      {
760        out =  Finley_RectangularMesh_Hex20_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info);
761            return out;
762      }    
763    
764      if( mpi_info->rank==0 )
765        domLeft = TRUE;
766      if( mpi_info->rank==mpi_info->size-1 )
767        domRight = TRUE;
768      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
769        domInternal = TRUE;
770    
771      /* dimensions of the local subdomain */
772      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );  
773    
774      NFaceElements=0;
775      if (!periodic[0]) {
776          NDOF0=N0;
777          NFaceElements+=(domRight+domLeft)*NE1*NE2;
778      } else {
779          NDOF0=N0-1;
780      }
781      if (!periodic[1]) {
782          NDOF1=N1;
783          NFaceElements+=2*numElementsLocal*NE2;
784      } else {
785          NDOF1=N1-1;
786      }
787      if (!periodic[2]) {
788          NDOF2=N2;
789          NFaceElements+=2*numElementsLocal*NE1;
790      } else {
791          NDOF2=N2-1;
792      }
793      
794        boundaryLeft = !domLeft || periodicLocal[0];
795        boundaryRight = !domRight || periodicLocal[1];
796        N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
797        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
798        firstNodeConstruct = firstNode - 2*boundaryLeft;
799        firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
800    
801      /*  allocate mesh: */
802      
803      sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
804      out=Finley_Mesh_alloc(name,3,order,mpi_info);
805    
806      if (! Finley_noError()) return NULL;
807    
808      out->Elements=Finley_ElementFile_alloc(Hex20,out->order,mpi_info);
809      if (useElementsOnFace) {
810         out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,mpi_info);
811         out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,mpi_info);
812      } else {
813         out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
814         out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,mpi_info);
815      }
816      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
817      if (! Finley_noError()) {
818          Finley_Mesh_dealloc(out);
819          return NULL;
820      }
821    
822      /*  allocate tables: */
823      Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
824      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*3, 0 );
825      Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1*NE2);
826      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
827      if (! Finley_noError()) {
828          Finley_Mesh_dealloc(out);
829          return NULL;
830      }
831    
832        k=0;
833      #pragma omp parallel for private(i0,i1,i2,k)
834      for (i2=0;i2<N2;i2++) {
835        for (i1=0;i1<N1;i1++) {
836          for (i0=0;i0<N0t;i0++,k++) {        
837            out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
838            out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
839            out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
840            out->Nodes->Id[k]=k;
841            out->Nodes->Tag[k]=0;
842            out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
843                    out->Nodes->Dom[k]=NODE_INTERNAL;
844          }
845        }
846      }
847    
848        /* mark the nodes that reference external and boundary DOF as such */
849        if( boundaryLeft ){
850            for( i1=0; i1<N1; i1++ )
851                for( i2=0; i2<N2; i2++ ) {
852                    out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
853                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_EXTERNAL;
854                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+2] = NODE_BOUNDARY;
855                }
856        }
857        if( boundaryRight ){
858            for( i1=0; i1<N1; i1++ )
859                for( i2=0; i2<N2; i2++ ) {
860                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
861                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
862                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-3] = NODE_BOUNDARY;
863                }
864        }
865        if( periodicLocal[0] ){
866            for( i1=0; i1<N1; i1++ )
867                for( i2=0; i2<N2; i2++ ) {
868                    out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
869                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
870                }
871        }
872            
873      /* tag Nodes that are referenced by face elements */
874      if (!periodic[2]) {    
875        for (i1=0;i1<N1;i1++) {
876          for (i0=0;i0<N0t;i0++) {  
877             out->Nodes->Tag[i0 + N0t*i1]+=100;
878             out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
879           }
880         }
881      }
882      if (!periodic[1]) {
883        for (i2=0;i2<N2;i2++) {
884          for (i0=0;i0<N0t;i0++) {
885             out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
886             out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
887          }
888        }
889      }
890      if (!periodic[0] && !domInternal ) {
891        for (i2=0;i2<N2;i2++) {
892          for (i1=0;i1<N1;i1++) {
893            if( domLeft )
894              out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
895            if( domRight )
896              out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
897          }
898        }
899      }
900    
901        /* form the boudary communication information */
902        forwardDOF  = MEMALLOC(NDOF1*NDOF2*2,index_t);
903        backwardDOF = MEMALLOC(NDOF1*NDOF2*2,index_t);
904        if( !(mpi_info->size==2 && periodicLocal[0])){
905            if( boundaryLeft  ) {
906                targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
907                for( iI=0, i2=0; i2<NDOF2; i2++ ){
908                    for( i1=0; i1<NDOF1; i1++, iI+=2 ){
909                        forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
910                        backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
911                        backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
912                    }
913                }
914                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
915                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
916            }
917            if( boundaryRight ) {
918                targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
919                for( iI=0,i2=0; i2<NDOF2; i2++ ){
920                    for( i1=0; i1<NDOF1; i1++, iI+=2 ){
921                        forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
922                        forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
923                        backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
924                    }
925                }
926                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
927                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
928            }
929        } else{
930            /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
931            targetDomain = 1;
932            
933            for( iI=0,i2=0; i2<NDOF2; i2++ ){
934                for( i1=0; i1<NDOF1; i1++, iI+=2 ){
935                    forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
936                    forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
937                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
938                }
939            }
940            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
941            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
942    
943            for( iI=0, i2=0; i2<NDOF2; i2++ ){
944                for( i1=0; i1<NDOF1; i1++, iI+=2 ){
945                    forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
946                    backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
947                    backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
948                }
949            }
950            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
951            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
952        }
953        MEMFREE( forwardDOF );
954        MEMFREE( backwardDOF );
955      /*   set the elements: */
956      
957        k=0;
958      #pragma omp parallel for private(i0,i1,i2,k,node0)
959      for (i2=0;i2<NE2;i2++) {
960        for (i1=0;i1<NE1;i1++) {
961          for (i0=0;i0<numElementsLocal;i0++,k++) {
962                    node0 = (periodicLocal[0] && !i0) ? 2*(i1*N0t + i2*N1*N0t) : 2*(i1*N0t + i2*N1*N0t + i0) + periodicLocal[0];
963    
964                    out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
965            out->Elements->Tag[k]=0;
966            out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
967            out->Elements->Dom[k]=ELEMENT_INTERNAL;
968    
969            out->Elements->Nodes[INDEX2(0,k,20)]=node0;
970            out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
971            out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0t+2;
972            out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0t;
973            out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0t*N1;
974            out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0t*N1+2;
975            out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0t*N1+2*N0t+2;
976            out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0t*N1+2*N0t;
977            out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
978            out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0t+2;
979            out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0t+1;
980            out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0t;
981            out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0t*N1;
982            out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0t*N1+2;
983            out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0t*N1+2*N0t+2;
984            out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0t*N1+2*N0t;
985            out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0t*N1+1;
986            out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0t*N1+N0t+2;
987            out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0t*N1+2*N0t+1;
988            out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0t*N1+N0t;
989          }
990        }
991      }
992      out->Elements->minColor=0;
993      out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
994    
995        if( boundaryLeft )
996            for( i2=0; i2<NE2; i2++ )
997                for( i1=0; i1<NE1; i1++ )
998                    out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
999        if( boundaryRight )
1000            for( i2=0; i2<NE2; i2++ )
1001                for( i1=0; i1<NE1; i1++ )
1002                    out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1003    
1004        out->Elements->numElements = numElementsLocal*NE1*NE2;
1005        Finley_ElementFile_setDomainFlags( out->Elements );
1006      
1007      /*   face elements: */
1008      
1009      if  (useElementsOnFace) {
1010         NUMNODES=20;
1011      } else {
1012         NUMNODES=8;
1013      }
1014      totalNECount=out->Elements->numElements;
1015      faceNECount=0;
1016        idCount = totalNECount;
1017      
1018      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
1019        numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1020      if (!periodic[0] && !domInternal) {
1021         /* **  elements on boundary 001 (x1=0): */
1022         if( domLeft ){
1023                 #pragma omp parallel for private(i1,i2,k)
1024                 for (i2=0;i2<NE2;i2++) {
1025                     for (i1=0;i1<NE1;i1++) {
1026                         k=i1+NE1*i2+faceNECount;
1027                         kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
1028                         facePerm =!useElementsOnFace ? face0R :  face0;
1029            
1030                         out->FaceElements->Id[k]=idCount++;
1031                         out->FaceElements->Tag[k]=1;
1032                         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1033                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
1034    
1035                         for( j=0; j<NUMNODES; j++ )
1036                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1037                     }
1038                 }
1039                 totalNECount+=NE1*NE2;
1040                 faceNECount+=NE1*NE2;
1041         }
1042         /* **  elements on boundary 002 (x1=1): */
1043             if( domRight ) {
1044                 #pragma omp parallel for private(i1,i2,k)
1045                 for (i2=0;i2<NE2;i2++) {
1046                     for (i1=0;i1<NE1;i1++) {
1047                         k=i1+NE1*i2+faceNECount;
1048                         kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
1049                         facePerm =!useElementsOnFace ? face1R :  face1;
1050            
1051                         out->FaceElements->Id[k]=idCount++;
1052                         out->FaceElements->Tag[k]=2;
1053                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1054                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
1055    
1056                         for( j=0; j<NUMNODES; j++ )
1057                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1058                     }
1059                 }
1060           totalNECount+=NE1*NE2;
1061           faceNECount+=NE1*NE2;
1062         }
1063      }
1064      if (!periodic[1]) {
1065         /* **  elements on boundary 010 (x2=0): */
1066      
1067         #pragma omp parallel for private(i0,i2,k)
1068         for (i2=0;i2<NE2;i2++) {
1069           for (i0=0;i0<numElementsLocal;i0++) {
1070             k=i0+numElementsLocal*i2+faceNECount;
1071                     kk=i0+numElementsLocal*NE1*i2;
1072                     facePerm =!useElementsOnFace ? face2R :  face2;
1073      
1074             out->FaceElements->Id[k]=idCount++;
1075             out->FaceElements->Tag[k]=10;
1076             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1077             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
1078    
1079                     for( j=0; j<NUMNODES; j++ )
1080                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1081           }
1082         }
1083           if( boundaryLeft ){
1084                for( i2=0; i2<NE2; i2++ )
1085                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1086                if( periodicLocal[0] )
1087                    for( i2=0; i2<NE2; i2++ )
1088                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1089             }
1090           if( boundaryRight )
1091                for( i2=0; i2<NE2; i2++ )
1092                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1093         totalNECount+=numElementsLocal*NE2;
1094         faceNECount+=numElementsLocal*NE2;
1095      
1096         /* **  elements on boundary 020 (x2=1): */
1097      
1098         #pragma omp parallel for private(i0,i2,k)
1099         for (i2=0;i2<NE2;i2++) {
1100           for (i0=0;i0<numElementsLocal;i0++) {
1101             k=i0+numElementsLocal*i2+faceNECount;
1102                     kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
1103                     facePerm =!useElementsOnFace ? face3R :  face3;
1104      
1105             out->FaceElements->Tag[k]=20;
1106             out->FaceElements->Id[k]=idCount++;
1107             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1108             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1109    
1110                     for( j=0; j<NUMNODES; j++ )
1111                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1112           }
1113         }
1114           if( boundaryLeft ){
1115                for( i2=0; i2<NE2; i2++ )
1116                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1117                if( periodicLocal[0] )
1118                    for( i2=0; i2<NE2; i2++ )
1119                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1120             }
1121           if( boundaryRight )
1122                for( i2=0; i2<NE2; i2++ )
1123                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1124         totalNECount+=numElementsLocal*NE2;
1125         faceNECount+=numElementsLocal*NE2;
1126      }
1127      if (!periodic[2]) {
1128         /*  elements on boundary 100 (x3=0): */
1129      
1130         #pragma omp parallel for private(i0,i1,k)
1131         for (i1=0;i1<NE1;i1++) {
1132           for (i0=0; i0<numElementsLocal; i0++) {
1133             k=i0+numElementsLocal*i1+faceNECount;
1134                     kk=i0 + i1*numElementsLocal;
1135                     facePerm =!useElementsOnFace ? face4R :  face4;
1136      
1137             out->FaceElements->Id[k]=idCount++;
1138             out->FaceElements->Tag[k]=100;
1139             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1140             out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
1141    
1142                     for( j=0; j<NUMNODES; j++ )
1143                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1144           }
1145         }
1146           if( boundaryLeft ){
1147                for( i1=0; i1<NE1; i1++ )
1148                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1149                if( periodicLocal[0] )
1150                    for( i1=0; i1<NE1; i1++ )
1151                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1152             }
1153           if( boundaryRight )
1154                for( i1=0; i1<NE1; i1++ )
1155                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1156         totalNECount+=NE1*numElementsLocal;
1157         faceNECount+=NE1*numElementsLocal;
1158            
1159         /* **  elements on boundary 200 (x3=1) */
1160      
1161         #pragma omp parallel for private(i0,i1,k)
1162         for (i1=0;i1<NE1;i1++) {
1163           for (i0=0;i0<numElementsLocal;i0++) {
1164             k=i0+numElementsLocal*i1+faceNECount;
1165                     kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
1166                     facePerm = !useElementsOnFace ? face5R : face5;
1167    
1168             out->FaceElements->Id[k]=idCount++;
1169             out->FaceElements->Tag[k]=200;
1170             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1171             out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
1172    
1173                     for( j=0; j<NUMNODES; j++ )
1174                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1175           }
1176         }
1177           if( boundaryLeft ){
1178                for( i1=0; i1<NE1; i1++ )
1179                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1180                if( periodicLocal[0] )
1181                    for( i1=0; i1<NE1; i1++ )
1182                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1183             }
1184           if( boundaryRight )
1185                for( i1=0; i1<NE1; i1++ )
1186                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1187         totalNECount+=NE1*numElementsLocal;
1188         faceNECount+=NE1*numElementsLocal;
1189      }
1190        out->FaceElements->elementDistribution->numInternal = faceNECount;
1191        
1192      out->FaceElements->minColor=0;
1193      out->FaceElements->maxColor=23;
1194        out->FaceElements->numElements=faceNECount;
1195        Finley_ElementFile_setDomainFlags( out->FaceElements );
1196    
1197      /* setup distribution info for other elements */
1198        Finley_ElementFile_setDomainFlags( out->ContactElements );
1199        Finley_ElementFile_setDomainFlags( out->Points );
1200    
1201        /* reorder the degrees of freedom */
1202        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1203    
1204      /*   condense the nodes: */
1205      Finley_Mesh_resolveNodeIds(out);
1206      if( !Finley_MPI_noError(mpi_info) )
1207      {
1208        Paso_MPIInfo_dealloc( mpi_info );
1209        Finley_Mesh_dealloc(out);
1210        return NULL;
1211      }
1212    
1213      /* prepare mesh for further calculatuions:*/
1214      Finley_Mesh_prepare(out);
1215      if( !Finley_MPI_noError(mpi_info) )
1216      {
1217        Paso_MPIInfo_dealloc( mpi_info );
1218        Finley_Mesh_dealloc(out);
1219        return NULL;
1220      }
1221    
1222      /* free up memory */
1223      Paso_MPIInfo_dealloc( mpi_info );
1224    
1225        //print_mesh_statistics( out, FALSE );
1226    
1227      #ifdef Finley_TRACE
1228      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1229      #endif
1230    
1231      return out;
1232    }
1233    #endif
1234    
Legend:
Removed from v.155  
changed lines
  Added in v.787
  ViewVC Help
Powered by ViewVC 1.1.26