/[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 1062 by gross, Mon Mar 26 06:17:53 2007 UTC revision 1312 by ksteube, Mon Sep 24 06:18:44 2007 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /* $Id$ */
3   *          Copyright 2006 by ACcESS MNRF                   *  
4   *                                                          *  /*******************************************************
5   *              http://www.access.edu.au                    *   *
6   *       Primary Business: Queensland, Australia            *   *           Copyright 2003-2007 by ACceSS MNRF
7   *  Licensed under the Open Software License version 3.0    *   *       Copyright 2007 by University of Queensland
8   *     http://www.opensource.org/licenses/osl-3.0.php       *   *
9   *                                                          *   *                http://esscc.uq.edu.au
10   ************************************************************   *        Primary Business: Queensland, Australia
11  */   *  Licensed under the Open Software License version 3.0
12     *     http://www.opensource.org/licenses/osl-3.0.php
13     *
14     *******************************************************/
15    
16  /**************************************************************/  /**************************************************************/
17    
# Line 19  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
 /*   Author: gross@access.edu.au */  
 /*   Version: $Id$ */  
   
 /**************************************************************/  
   
25  #include "RectangularMesh.h"  #include "RectangularMesh.h"
26    
 #ifdef PASO_MPI  
 /* get the number of nodes/elements for domain with rank=rank, of size processors  
    where n is the total number of nodes/elements in the global domain */  
 static index_t domain_MODdim( index_t rank, index_t size, index_t n )  
 {  
   rank = size-rank-1;  
   
   if( rank < n%size )  
     return (index_t)floor(n/size)+1;  
   return (index_t)floor(n/size);  
 }  
   
27    
28  /* Determines the number of nodes/elements etc along an axis which is numElementsGlobal long for domain rank */  Finley_Mesh* Finley_RectangularMesh_Rec8(dim_t* numElements,
29  /* A bit messy, but it only has to be done once... */                                            double* Length,
30  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 )                                            bool_t* periodic,
31                                              index_t order,
32                                              index_t reduced_order,
33                                              bool_t useElementsOnFace,
34                                              bool_t useFullElementOrder,
35                                              bool_t optimize)
36  {  {
37    index_t i0;    #define N_PER_E 2
38    dim_t numNodesGlobal = numElementsGlobal+1;    #define DIM 2
39      dim_t N0,N1,NE0,NE1,i0,i1,k,Nstride0,Nstride1;
40      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NN, local_NE0, local_NE1, local_N0, local_N1;
41      index_t e_offset1, e_offset0, offset0, offset1, global_i0, global_i1;
42      index_t node0, myRank;
43      Finley_Mesh* out;
44      Paso_MPIInfo *mpi_info = NULL;
45      char name[50];
46      double time0=Finley_timer();
47    
48    (*numNodesLocal) = 2*domain_MODdim( rank, size, numNodesGlobal ) - 1;    /* get MPI information */
49    if( rank<size-1 ) // add on node for right hand boundary    mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
50      (*numNodesLocal) += 1;    if (! Finley_noError()) {
51            return NULL;
   numElementsLocal[0] = domain_MODdim( rank, size, numNodesGlobal )+1;  
   periodicLocal[0] = periodicLocal[1] = FALSE;  
   nodesExternal[0] = nodesExternal[1] = 1;  
   if( periodic )  
   {  
     if( size==1 )  
     {  
       numElementsLocal[0] = numElementsGlobal;  
       nodesExternal[0] = nodesExternal[1] = 0;  
       periodicLocal[0] = periodicLocal[1] = TRUE;  
     }  
     else  
     {  
       if( rank==0 )  
       {  
         periodicLocal[0] = TRUE;  
         numNodesLocal[0]++;  
       }  
       if( rank==(size-1) )  
       {        
         periodicLocal[1] = TRUE;  
         numNodesLocal[0]--;  
         numElementsLocal[0]--;  
       }  
     }  
   }  
   else if( !periodic )  
   {  
     if( rank==0 ){  
       nodesExternal[0]--;  
       numElementsLocal[0]--;  
     }  
     if( rank==(size-1) )  
     {  
       nodesExternal[1]--;  
       numElementsLocal[0]--;  
     }  
   }  
   nodesExternal[0]*=2;  
   numNodesExternal[0] = nodesExternal[0]+nodesExternal[1];  
     
   numElementsInternal[0] = numElementsLocal[0];  
   if( (rank==0) && (rank==size-1) );  
       
   else if( !periodic && ( (rank==0) ^ (rank==size-1) ) )  
       numElementsInternal[0] -= 1;  
   else  
     numElementsInternal[0] -= 2;  
   
   firstNode[0] = 0;  
   for( i0=0; i0<rank; i0++ )  
     firstNode[0] += domain_MODdim( i0, size, numNodesGlobal );  
   firstNode[0] *= 2;  
     
   numDOFLocal[0] = numNodesLocal[0];  
   if( periodicLocal[0] )  
     numDOFLocal[0]--;  
     
   DOFExternal[0] = nodesExternal[0];  
   DOFExternal[1] = nodesExternal[1];  
   
   if( size>1 )  
   {  
     DOFBoundary[0] = periodicLocal[0]*1 + (rank>0)*1;  
     DOFBoundary[1] = periodicLocal[1]*2 + (rank<(size-1))*2;  
   }  
   else  
   {  
     DOFBoundary[0] = DOFBoundary[1] = 0;  
52    }    }
53  }    myRank=mpi_info->rank;
   
54    
55  #endif    /* set up the global dimensions of the mesh */
 /**************************************************************/  
56    
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Rec8_singleCPU(int* numElements,double* Length,int* periodic,int order, index_t reduced_order, int useElementsOnFace, Paso_MPIInfo *mpi_info)  
 #else  
 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order, index_t reduced_order, int useElementsOnFace)  
 #endif  
 {  
   dim_t N0,N1,NE0,NE1,i0,i1,totalNECount,faceNECount,NDOF0,NDOF1,NFaceElements,NUMNODES,M0,M1;  
   index_t k,node0;  
   Finley_Mesh* out=NULL;  
   char name[50];  
   double time0=Finley_timer();  
57    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
58    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
59    N0=2*NE0+1;    N0=N_PER_E*NE0+1;
60    N1=2*NE1+1;    N1=N_PER_E*NE1+1;
   
 #ifndef PASO_MPI  
   if (N0<=N1) {  
      M0=1;  
      M1=N0;  
   } else {  
      M0=N1;  
      M1=1;  
   }  
 #else  
      M0=1;  
      M1=N0;  
 #endif  
61    
62    NFaceElements=0;    /*  allocate mesh: */  
   if (!periodic[0]) {  
       NDOF0=N0;  
       NFaceElements+=2*NE1;  
   } else {  
       NDOF0=N0-1;  
   }  
   if (!periodic[1]) {  
       NDOF1=N1;  
       NFaceElements+=2*NE0;  
   } else {  
       NDOF1=N1-1;  
   }  
   
   /*  allocate mesh: */  
     
63    sprintf(name,"Rectangular %d x %d mesh",N0,N1);    sprintf(name,"Rectangular %d x %d mesh",N0,N1);
64  #ifdef PASO_MPI    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
65    out=Finley_Mesh_alloc(name,2,order, reduced_order,mpi_info);    if (! Finley_noError()) {
66  #else        Paso_MPIInfo_free( mpi_info );
   out=Finley_Mesh_alloc(name,2,order, reduced_order);  
 #endif  
   if (! Finley_noError()) return NULL;  
   
 #ifdef PASO_MPI  
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order,mpi_info);  
 #else  
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order);  
 #endif  
   
   if (! Finley_noError()) {  
        Finley_Mesh_dealloc(out);  
        return NULL;  
   }  
     
   /*  allocate tables: */  
   Finley_NodeFile_allocTable(out->Nodes,N0*N1);  
   Finley_ElementFile_allocTable(out->Elements,NE0*NE1);  
   Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);  
 #ifdef PASO_MPI  
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1, 0, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1, NE0*NE1);  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );  
 #endif  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
67        return NULL;        return NULL;
68    }    }
69    
70    /*  set nodes: */    if (useFullElementOrder) {
71    #pragma omp parallel for private(i0,i1,k)       Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec9,
72    for (i1=0;i1<N1;i1++) {                                              out->order,
73      for (i0=0;i0<N0;i0++) {                                              out->reduced_order,
74        k=M0*i0+M1*i1;                                              mpi_info));
75        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE(i0)/DBLE(N0-1)*Length[0];       if (useElementsOnFace) {
76        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];           Finley_setError(SYSTEM_ERROR,"rich elements for Rec9 elements is not supported yet.");
77        out->Nodes->Id[k]=i0+N0*i1;       } else {
78        out->Nodes->Tag[k]=0;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
79        out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1);                                                      out->order,
80  #ifdef PASO_MPI                                                      out->reduced_order,
81        out->Nodes->Dom[k]=NODE_INTERNAL;                                                      mpi_info));
82  #endif           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
83      }                                                         out->order,
84    }                                                         out->reduced_order,
85    /* tags for the faces: */                                                         mpi_info));
   if (!periodic[1]) {  
     for (i0=0;i0<N0;i0++) {  
       out->Nodes->Tag[M0*i0+M1*0]+=10;  
       out->Nodes->Tag[M0*i0+M1*(N1-1)]+=20;  
     }  
   }  
   if (!periodic[0]) {  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[M0*0+M1*i1]+=1;  
       out->Nodes->Tag[M0*(N0-1)+M1*i1]+=2;  
     }  
   }  
   
   /*   set the elements: */  
     
   #pragma omp parallel for private(i0,i1,k,node0)  
   for (i1=0;i1<NE1;i1++) {  
     for (i0=0;i0<NE0;i0++) {  
       k=i0+NE0*i1;  
       node0=2*i0+2*i1*N0;  
   
       out->Elements->Id[k]=k;  
       out->Elements->Tag[k]=0;  
       out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);  
 #ifdef PASO_MPI  
       out->Elements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
       out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;  
       out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0+2;  
       out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0;  
       out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;  
       out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0+2;  
       out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0+1;  
       out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0;  
   
     }  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     
   /*   face elements: */  
   if (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=3;  
   }  
     
   totalNECount=NE0*NE1;  
   faceNECount=0;  
     
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x2=0): */  
     
      #pragma omp parallel for private(i0,k,node0)  
      for (i0=0;i0<NE0;i0++) {  
        k=i0+faceNECount;  
        node0=2*i0;  
   
        out->FaceElements->Id[k]=i0+totalNECount;  
        out->FaceElements->Tag[k]=10;  
        out->FaceElements->Color[k]=i0%2;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+1;  
        }  
86       }       }
87       totalNECount+=NE0;    } else  {
88       faceNECount+=NE0;       Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Rec8,out->order,out->reduced_order,mpi_info));
89           if (useElementsOnFace) {
90       /* **  elements on boundary 020 (x2=1): */           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8Face,
91                                                                        out->order,
92       #pragma omp parallel for private(i0,k,node0)                                                                    out->reduced_order,
93       for (i0=0;i0<NE0;i0++) {                                                                    mpi_info));
94         k=i0+faceNECount;           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8Face_Contact,
95         node0=2*i0+2*(NE1-1)*N0;                                                                      out->order,
96                                                                        out->reduced_order,
97         out->FaceElements->Id[k]=i0+totalNECount;                                                                      mpi_info));
98         out->FaceElements->Tag[k]=20;       } else {
99         out->FaceElements->Color[k]=i0%2+2;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Line3,
100  #ifdef PASO_MPI                                                                    out->order,
101           out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                                                                    out->reduced_order,
102  #endif                                                                    mpi_info));
103         if (useElementsOnFace) {           Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Line3_Contact,
104            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;                                                                       out->order,
105            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;                                                                       out->reduced_order,
106            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0;                                                                       mpi_info));
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+2;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+1;  
        }  
107       }       }
      totalNECount+=NE0;  
      faceNECount+=NE0;  
108    }    }
109    if (!periodic[0]) {    Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
110       /* **  elements on boundary 001 (x1=0): */                                                   out->order,
111                                                       out->reduced_order,
112       #pragma omp parallel for private(i1,k,node0)                                                   mpi_info));
      for (i1=0;i1<NE1;i1++) {  
        k=i1+faceNECount;  
        node0=2*i1*N0;  
   
        out->FaceElements->Id[k]=i1+totalNECount;  
        out->FaceElements->Tag[k]=1;  
        out->FaceElements->Color[k]=(i1%2)+4;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
   
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0;  
        }  
     
      }  
      totalNECount+=NE1;  
      faceNECount+=NE1;  
     
      /* **  elements on boundary 002 (x1=1): */  
       
      #pragma omp parallel for private(i1,k,node0)  
      for (i1=0;i1<NE1;i1++) {  
        k=i1+faceNECount;  
        node0=2*(NE0-1)+2*i1*N0;  
     
        out->FaceElements->Id[k]=i1+totalNECount;  
        out->FaceElements->Tag[k]=2;  
        out->FaceElements->Color[k]=(i1%2)+6;  
 #ifdef PASO_MPI  
          out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
 #endif  
     
        if (useElementsOnFace) {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0;  
           out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0;  
           out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;  
           out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;  
           out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0;  
           out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;  
        } else {  
           out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;  
           out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;  
           out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+2;  
        }  
      }  
      totalNECount+=NE1;  
      faceNECount+=NE1;  
   
   }  
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=7;  
   
 #ifdef PASO_MPI  
     Finley_ElementFile_setDomainFlags( out->Elements );  
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     /* reorder the degrees of freedom */  
     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );  
 #endif  
   
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
       
113    if (! Finley_noError()) {    if (! Finley_noError()) {
114        Finley_Mesh_dealloc(out);        Paso_MPIInfo_free( mpi_info );
115          Finley_Mesh_free(out);
116        return NULL;        return NULL;
117    }    }
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
118    
119      /* work out the largest dimension */
120      if (N1==MAX(N0,N1)) {
121         Nstride0=1;
122         Nstride1=N0;
123         local_NE0=NE0;
124         e_offset0=0;
125         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
126      } else {
127         Nstride0=N1;
128         Nstride1=1;
129         Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
130         local_NE1=NE1;
131         e_offset1=0;
132      }
133      offset0=e_offset0*N_PER_E;
134      offset1=e_offset1*N_PER_E;
135      local_N0=local_NE0*N_PER_E+1;
136      local_N1=local_NE1*N_PER_E+1;
137    
138    #ifdef Finley_TRACE    /* get the number of surface elements */
   printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);  
   #endif  
   
   if (Finley_noError()) {  
        if ( ! Finley_Mesh_isPrepared(out)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   } else {  
       Finley_Mesh_dealloc(out);  
   }  
   return out;  
 }  
   
 #ifdef PASO_MPI  
 Finley_Mesh* Finley_RectangularMesh_Rec8(int* numElements,double* Length,int* periodic,int order, index_t reduced_order ,int useElementsOnFace) {  
   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];  
     dim_t N0t, NDOF0t;  
   index_t *numForward=NULL, *numBackward=NULL;  
   index_t kk, k,i,j,node0,firstNodeConstruct,firstNode=0, DOFcount=0, forwardDOF=NULL, backwardDOF=NULL, targetDomain=-1,faceNEBoundary;  
   Finley_Mesh* out=NULL;  
   char name[50];  
     index_t *indexForward=NULL, *indexBackward=NULL;  
   double time0=Finley_timer();  
   bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;  
   Paso_MPIInfo *mpi_info = NULL;  
     index_t faceStart, numFaceLeft;  
   
     /* these are used in constructing the face elements, and give the permutation  
        of the node order of a normal element for each face */  
     index_t *facePerm;    
     index_t face0[]={0, 1, 2, 3, 4, 5, 6, 7};  
     index_t face1[]={2, 3, 0, 1, 6, 7, 4, 5};  
     index_t face2[]={3, 0, 1, 2, 7, 4, 5, 6};  
     index_t face3[]={1, 2, 3, 0, 5, 6, 7, 4};  
     index_t face0nc[]={0, 1, 4};  
     index_t face1nc[]={2, 3, 6};  
     index_t face2nc[]={3, 0, 7};  
     index_t face3nc[]={1, 2, 5};  
   
   NE0=MAX(1,numElements[0]);  
   NE1=MAX(1,numElements[1]);  
   N0=2*NE0+1;  
   N1=2*NE1+1;  
   
   /* get MPI information */  
   mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );  
   if (! Finley_noError()) {  
         Finley_Mesh_dealloc(out);  
         return NULL;  
   }  
   
     // use the serial method for the single CPU case  
     if( mpi_info->size==1 ){  
         out = Finley_RectangularMesh_Rec8_singleCPU( numElements,Length,periodic,order, reduced_order,useElementsOnFace, mpi_info);  
         return out;  
     }  
       
   if( mpi_info->rank==0 )  
     domLeft = TRUE;  
   if( mpi_info->rank==mpi_info->size-1 )  
     domRight = TRUE;  
   if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )  
     domInternal = TRUE;  
   
   /* dimensions of the local subdomain */  
   domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );    
139    
140    NFaceElements=0;    NFaceElements=0;
141    if (!periodic[0]) {    if (!periodic[0]) {
142        NDOF0=N0;       NDOF0=N0;
143        NFaceElements+=(domLeft + domRight)*NE1;       if (e_offset0 == 0) NFaceElements+=local_NE1;
144         if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1;
145    } else {    } else {
146        NDOF0=N0-1;        NDOF0=N0-1;
147    }    }
148    if (!periodic[1]) {    if (!periodic[1]) {
149        NDOF1=N1;       NDOF1=N1;
150        NFaceElements+=2*numElementsLocal;       if (e_offset1 == 0) NFaceElements+=local_NE0;
151         if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0;
152    } else {    } else {
153        NDOF1=N1-1;        NDOF1=N1-1;
154    }    }
155    
     boundaryLeft = !domLeft || periodicLocal[0];  
     boundaryRight = !domRight || periodicLocal[1];  
     N0t = numNodesLocal + boundaryRight + boundaryLeft*2;  
     NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;  
     firstNodeConstruct = firstNode - 2*boundaryLeft;  
     firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;  
   
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d mesh",N0,N1);  
   out=Finley_Mesh_alloc(name,2,order,reduced_order,mpi_info);  
   
   if (! Finley_noError()) return NULL;  
   
   out->Elements=Finley_ElementFile_alloc(Rec8,out->order, out->reduced_order,mpi_info);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Rec8Face,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Rec8Face_Contact,out->order, out->reduced_order,mpi_info);  
   } else {  
      out->FaceElements=Finley_ElementFile_alloc(Line3,out->order, out->reduced_order,mpi_info);  
      out->ContactElements=Finley_ElementFile_alloc(Line3_Contact,out->order, out->reduced_order,mpi_info);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order, out->reduced_order,mpi_info);  
   
   if (! Finley_noError()) {  
        Finley_Mesh_dealloc(out);  
        return NULL;  
   }  
     
156    /*  allocate tables: */    /*  allocate tables: */
157    Finley_NodeFile_allocTable(out->Nodes,N0t*N1);  
158    Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1);    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1);
159      Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1);
160    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
       
   Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0t*NDOF1, NDOF1*3, 0 );  
   Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1, NE1*(numElementsLocal-boundaryRight*(!periodic[1])) );  
   Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(!periodic[1]) );  
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
161    
162    /*  set nodes: */    if (Finley_noError()) {
163    #pragma omp parallel for private(i0,i1,k)       /* create nodes */
164    for (i1=0;i1<N1;i1++) {    
165      for (i0=0;i0<N0t;i0++, k++) {       #pragma omp parallel for private(i0,i1,k,global_i0,global_i1)
166              k=i0+i1*N0t;       for (i1=0;i1<local_N1;i1++) {
167        out->Nodes->Coordinates[INDEX2(0,k,2)]=DBLE((firstNodeConstruct + i0) % N0)/DBLE(N0-1)*Length[0];         for (i0=0;i0<local_N0;i0++) {
168        out->Nodes->Coordinates[INDEX2(1,k,2)]=DBLE(i1)/DBLE(N1-1)*Length[1];             k=i0+local_N0*i1;
169        out->Nodes->Id[k]=i0 + i1*N0t;             global_i0=i0+offset0;
170        out->Nodes->Tag[k]=0;             global_i1=i1+offset1;
171        out->Nodes->Dom[k]=NODE_INTERNAL;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
172        out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
173      }             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1;
174          /* take care to make the DOF reflect the internal/external position of the DOF */             out->Nodes->Tag[k]=0;
175          if( boundaryLeft ){             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
176              out->Nodes->Dom[i1*N0t] = NODE_EXTERNAL;                                                     +Nstride1*(global_i1%NDOF1);
177              out->Nodes->Dom[i1*N0t+1] = NODE_EXTERNAL;           }
             out->Nodes->Dom[i1*N0t+2] = NODE_BOUNDARY;    
             if( periodicLocal[0] ){  
                 out->Nodes->degreeOfFreedom[i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];  
                 out->Nodes->Dom[i1*N0t+3] = NODE_BOUNDARY;    
             }  
         }  
         if( boundaryRight ){  
             out->Nodes->Dom[(i1+1)*N0t-1] = NODE_EXTERNAL;    
             out->Nodes->Dom[(i1+1)*N0t-2] = NODE_BOUNDARY;    
             out->Nodes->Dom[(i1+1)*N0t-3] = NODE_BOUNDARY;    
         }  
   }  
   /* tags for the faces: */  
   if (!periodic[1]) {  
     for (i0=0;i0<N0t;i0++) {  
       out->Nodes->Tag[i0]+=10;  
       out->Nodes->Tag[i0+N0t*(N1-1)]+=20;  
     }  
   }  
   if (domLeft && !periodic[0]) {  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[i1*N0t]+=1;  
     }  
   }  
     if (domRight && !periodic[0]){  
     for (i1=0;i1<N1;i1++) {  
       out->Nodes->Tag[(i1+1)*N0t-1]+=2;  
         }  
     }  
     
     /* setup the receive information */  
     indexBackward = TMPMEMALLOC( NDOF1*2, index_t );  
     indexForward  = TMPMEMALLOC( NDOF1*2, index_t );  
     if( !(mpi_info->size==2 && periodicLocal[0])){  
         /* Backward (receive) DOF */  
         if( boundaryLeft ){  
             targetDomain = mpi_info->rank-1<0 ? mpi_info->size-1 : mpi_info->rank-1;  
   
             for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];  
                 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];  
                 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];  
             }  
   
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );  
         }  
           
         if( boundaryRight ){  
             targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;  
   
             for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];  
                 indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];  
                 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];  
             }  
   
             Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );  
       Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );  
         }  
     }  
     else{  
         targetDomain = 1;  
   
         for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i] = out->Nodes->degreeOfFreedom[(i+1)*N0t-1];  
                 indexForward[i1]   = out->Nodes->degreeOfFreedom[(i+1)*N0t-3];  
                 indexForward[i1+1] = out->Nodes->degreeOfFreedom[(i+1)*N0t-2];  
         }  
   
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexForward );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexBackward );  
   
         for( i=0, i1=0; i<NDOF1; i++, i1+=2 ){  
                 indexBackward[i1] = out->Nodes->degreeOfFreedom[i*N0t];  
                 indexBackward[i1+1] = out->Nodes->degreeOfFreedom[i*N0t+1];  
                 indexForward[i] = out->Nodes->degreeOfFreedom[i*N0t+2];  
         }  
   
         Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1, indexForward );  
         Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*2, indexBackward );  
     }  
     TMPMEMFREE( indexBackward );  
     TMPMEMFREE( indexForward );  
   
   /*   set the elements: */  
   #pragma omp parallel for private(i0,i1,k,node0)  
   for (i1=0;i1<NE1;i1++) {  
     for (i0=0;i0<numElementsLocal;i0++) {  
       k=i0+numElementsLocal*i1;  
             node0 = (periodicLocal[0] && !i0) ? 2*i1*N0t :  2*i0+2*i1*N0t+periodicLocal[0] ;  
   
       out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0) * NE1 + i1;  
       out->Elements->Tag[k]=0;  
       out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1);  
             out->Elements->Dom[k] = ELEMENT_INTERNAL;  
   
       out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
       out->Elements->Nodes[INDEX2(1,k,8)]=node0+2;  
       out->Elements->Nodes[INDEX2(2,k,8)]=node0+2*N0t+2;  
       out->Elements->Nodes[INDEX2(3,k,8)]=node0+2*N0t;  
       out->Elements->Nodes[INDEX2(4,k,8)]=node0+1;  
       out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0t+2;  
       out->Elements->Nodes[INDEX2(6,k,8)]=node0+2*N0t+1;  
       out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0t;  
     }  
         if( boundaryLeft )  
             out->Elements->Dom[i1*numElementsLocal] = ELEMENT_BOUNDARY;  
         if( boundaryRight )  
             out->Elements->Dom[(i1+1)*numElementsLocal-1] = ELEMENT_BOUNDARY;  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0);  
     
     /* set the distribution information for the elements */  
     Finley_ElementFile_setDomainFlags( out->Elements );  
   
   /*   face elements: */  
   if (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=3;  
   }  
     
     Finley_ElementFile_setDomainFlags( out->FaceElements );  
   totalNECount=numElementsLocal*NE1;  
   faceNECount=0;  
     faceNEBoundary = totalNECount + numElementsInternal*(!periodic[1])*2 + NE1*(!periodic[0])*(domLeft+domRight);  
     
     numFaceLeft = domLeft*(!periodic[0])*NE1;  
     faceStart = out->FaceElements->elementDistribution->vtxdist[mpi_info->rank];  
   if (!periodic[1]) {  
      /* **  elements on boundary 010 (x1=0): */  
     
      #pragma omp parallel for private(i0,k,kk,facePerm)  
      for (i0=0;i0<numElementsLocal; i0++) {  
        k=i0+faceNECount;  
              kk=i0;  
              facePerm = useElementsOnFace ? face0 : face0nc;  
   
        out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2;  
        out->FaceElements->Tag[k]=10;  
        out->FaceElements->Color[k]=i0%2;  
              out->FaceElements->Dom[k]=ELEMENT_INTERNAL;  
                   
              for( j=0; j<NUMNODES; j++ )  
                   out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
178       }       }
179           if( boundaryLeft ){       /*   set the elements: */
180               out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;       NN=out->Elements->numNodes;
181               if( periodicLocal[0] )       #pragma omp parallel for private(i0,i1,k,node0)
182                   out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;       for (i1=0;i1<local_NE1;i1++) {
183           }           for (i0=0;i0<local_NE0;i0++) {
184           if( boundaryRight ){            
185               out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;             k=i0+local_NE0*i1;        
186               out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1];             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
187           }  
188       totalNECount+=numElementsLocal;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1);
189       faceNECount+=numElementsLocal;             out->Elements->Tag[k]=0;
190                 out->Elements->Owner[k]=myRank;
191       /* **  elements on boundary 020 (x1=1): */  
192                 out->Elements->Nodes[INDEX2(0,k,NN)]=node0;
193       #pragma omp parallel for private(i0,k,kk,facePerm)             out->Elements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
194       for (i0=0;i0<numElementsLocal;i0++) {             out->Elements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
195         k=i0+faceNECount;             out->Elements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
196               kk=i0 + numElementsLocal*(NE1-1);             out->Elements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
197               facePerm = useElementsOnFace ? face1 : face1nc;             out->Elements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
198               out->Elements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
199         out->FaceElements->Id[k]=faceStart + numFaceLeft + i0*2+1;             out->Elements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
200         out->FaceElements->Tag[k]=20;             if (useFullElementOrder) {
201         out->FaceElements->Color[k]=i0%2+2;                out->Elements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride1+1*Nstride0;
202               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;             }
203             }
              for( j=0; j<NUMNODES; j++ )  
                   out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];  
204       }       }
205           if( boundaryLeft ){       /* face elements */
206               out->FaceElements->Dom[faceNECount]=ELEMENT_BOUNDARY;       NN=out->FaceElements->numNodes;
207               if( periodicLocal[0] )       totalNECount=NE0*NE1;
208                   out->FaceElements->Dom[faceNECount+1]=ELEMENT_BOUNDARY;       faceNECount=0;
209           }       if (!periodic[0]) {
210           if( boundaryRight ){          /* **  elements on boundary 001 (x1=0): */
211               out->FaceElements->Dom[faceNECount+numElementsLocal-1]=ELEMENT_BOUNDARY;      
212               out->FaceElements->Id[faceNECount+numElementsLocal-1]=out->FaceElements->elementDistribution->vtxdist[mpi_info->rank+1]+1;          if (e_offset0 == 0) {
213           }             #pragma omp parallel for private(i1,k,node0)
214       totalNECount+=numElementsLocal;             for (i1=0;i1<local_NE1;i1++) {
215       faceNECount+=numElementsLocal;        
216    }                 k=i1+faceNECount;
217    if (domLeft && !periodicLocal[0]) {                 node0=Nstride1*N_PER_E*(i1+e_offset1);
218       /* **  elements on boundary 001 (x0=0): */  
219                     out->FaceElements->Id[k]=i1+e_offset1+totalNECount;
220       #pragma omp parallel for private(i1,k,kk,facePerm)                 out->FaceElements->Tag[k]=1;
221       for (i1=0;i1<NE1;i1++) {                 out->FaceElements->Owner[k]=myRank;
222         k=i1+faceNECount;                 if (useElementsOnFace) {
223               kk=i1*numElementsLocal;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
224               facePerm = useElementsOnFace ? face2 : face2nc;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
225                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride0;
226         out->FaceElements->Id[k]=faceStart + i1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0;
227         out->FaceElements->Tag[k]=1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1;
228         out->FaceElements->Color[k]=(i1%2)+4;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride0;
229               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1+2*Nstride0;
230                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride1+1*Nstride0;
231               for( j=0; j<NUMNODES; j++ )                 } else {
232                    out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1;
233                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0;
234                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1;
235                   }
236               }
237               faceNECount+=local_NE1;
238            }
239            totalNECount+=NE1;
240            /* **  elements on boundary 002 (x1=1): */
241            if (local_NE0+e_offset0 == NE0) {
242               #pragma omp parallel for private(i1,k,node0)
243               for (i1=0;i1<local_NE1;i1++) {
244                   k=i1+faceNECount;
245                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1);
246    
247                   out->FaceElements->Id[k]=(i1+e_offset1)+totalNECount;
248                   out->FaceElements->Tag[k]=2;
249                   out->FaceElements->Owner[k]=myRank;
250    
251                   if (useElementsOnFace) {
252                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
253                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
254                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1;
255                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0;
256                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+Nstride1+2*Nstride0;
257                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+1*Nstride0;
258                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride1;
259                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride0;
260                   } else {
261                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride0;
262                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1+2*Nstride0;
263                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride1+2*Nstride0;
264                   }
265               }
266               faceNECount+=local_NE1;
267             }
268             totalNECount+=NE1;
269       }       }
270       totalNECount+=NE1;       if (!periodic[1]) {
271       faceNECount+=NE1;          /* **  elements on boundary 010 (x2=0): */
272              if (e_offset1 == 0) {
273       /* **  elements on boundary 002 (x0=1): */             #pragma omp parallel for private(i0,k,node0)
274    }             for (i0=0;i0<local_NE0;i0++) {
275    if (domRight && !periodicLocal[1]) {                 k=i0+faceNECount;
276       #pragma omp parallel for private(i1,k,kk,facePerm)                 node0=Nstride0*N_PER_E*(i0+e_offset0);
277       for (i1=0;i1<NE1;i1++) {          
278         k=i1+faceNECount;                 out->FaceElements->Id[k]=e_offset0+i0+totalNECount;
279               kk=(i1+1)*numElementsLocal-1;                 out->FaceElements->Tag[k]=10;
280               facePerm = useElementsOnFace ? face3 : face3nc;                 out->FaceElements->Owner[k]=myRank;
281          
282                   if (useElementsOnFace) {
283                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
284                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
285                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+2*Nstride0;
286                       out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1;
287                       out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride0;
288                       out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+2*Nstride0;
289                       out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride1+1*Nstride0;
290                       out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1;
291                   } else {
292                       out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0;
293                       out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride0;
294                       out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+1*Nstride0;
295                   }
296               }
297               faceNECount+=local_NE0;
298            }
299            totalNECount+=NE0;
300            /* **  elements on boundary 020 (x2=1): */
301            if (local_NE1+e_offset1 == NE1) {
302               #pragma omp parallel for private(i0,k,node0)
303               for (i0=0;i0<local_NE0;i0++) {
304                   k=i0+faceNECount;
305                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1);
306        
307         out->FaceElements->Id[k]=faceStart + numElementsLocal*2*(!periodic[1])+i1;                 out->FaceElements->Id[k]=i0+e_offset0+totalNECount;
308         out->FaceElements->Tag[k]=2;                 out->FaceElements->Tag[k]=20;
309         out->FaceElements->Color[k]=(i1%2)+6;                 out->FaceElements->Owner[k]=myRank;
310               out->FaceElements->Dom[k]=ELEMENT_INTERNAL;                 if (useElementsOnFace) {
311                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
312                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
313                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0;
314                        out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride0;
315                        out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+2*Nstride1+1*Nstride0;
316                        out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1;
317                        out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride0;
318                        out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride1+2*Nstride0;
319                   } else {
320                        out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+2*Nstride1+2*Nstride0;
321                        out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride1;
322                        out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride1+1*Nstride0;
323                   }
324               }
325               faceNECount+=local_NE0;
326            }
327            totalNECount+=NE0;
328         }
329         /* add tag names */
330         Finley_Mesh_addTagMap(out,"top", 20);
331         Finley_Mesh_addTagMap(out,"bottom", 10);
332         Finley_Mesh_addTagMap(out,"left", 1);
333         Finley_Mesh_addTagMap(out,"right", 2);
334        
335               for( j=0; j<NUMNODES; j++ )       /* prepare mesh for further calculatuions:*/
336                    out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,8)];       if (Finley_noError()) {
337             Finley_Mesh_resolveNodeIds(out);
338         }
339         if (Finley_noError()) {
340             Finley_Mesh_prepare(out, optimize);
341       }       }
      totalNECount+=NE1;  
      faceNECount+=NE1;  
342    }    }
343    out->FaceElements->minColor=0;    if (!Finley_noError()) {
344    out->FaceElements->maxColor=7;        Finley_Mesh_free(out);
     out->FaceElements->numElements = faceNECount;  
   
   /* setup distribution info for other elements */  
     Finley_ElementFile_setDomainFlags( out->ContactElements );  
     Finley_ElementFile_setDomainFlags( out->Points );  
   
     Finley_Mesh_prepareElementDistribution( out );  
   
     /* reorder the degrees of freedom */  
     Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );  
       
   /*   condense the nodes: */  
   Finley_Mesh_resolveNodeIds(out);  
   if ( !Finley_MPI_noError( mpi_info )) {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
345    }    }
346      /* free up memory */
347    /* prepare mesh for further calculatuions:*/    Paso_MPIInfo_free( mpi_info );  
   Finley_Mesh_prepare(out) ;  
   
348    #ifdef Finley_TRACE    #ifdef Finley_TRACE
349    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
350    #endif    #endif
351    
   if ( !Finley_MPI_noError( mpi_info )) {  
     Paso_MPIInfo_dealloc( mpi_info );  
     Finley_Mesh_dealloc(out);  
     return NULL;  
   }  
   if (Finley_noError()) {  
        if (! Finley_Mesh_isPrepared(out)) {  
           Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");  
        }  
   }  
352    return out;    return out;
353  }  }
 #endif  

Legend:
Removed from v.1062  
changed lines
  Added in v.1312

  ViewVC Help
Powered by ViewVC 1.1.26