/[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/esys2/finley/src/finleyC/Mesh_hex20.c revision 150 by jgs, Thu Sep 15 03:44:45 2005 UTC trunk/finley/src/Mesh_hex20.c revision 964 by gross, Tue Feb 13 05:10:26 2007 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    dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;  #endif
133      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;
136    char name[50];    char name[50];
# Line 46  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 143  Finley_Mesh* Finley_RectangularMesh_Hex2
143    N1=2*NE1+1;    N1=2*NE1+1;
144    N2=2*NE2+1;    N2=2*NE2+1;
145    
146      if (N0<=MIN(N1,N2)) {
147         if (N1 <= N2) {
148            M0=1;
149            M1=N0;
150            M2=N0*N1;
151         } else {
152            M0=1;
153            M2=N0;
154            M1=N0*N2;
155         }
156      } else if (N1<=MIN(N2,N0)) {
157         if (N2 <= N0) {
158            M1=1;
159            M2=N1;
160            M0=N2*N1;
161         } else {
162            M1=1;
163            M0=N1;
164            M2=N1*N0;
165         }
166      } else {
167         if (N0 <= N1) {
168            M2=1;
169            M0=N2;
170            M1=N2*N0;
171         } else {
172            M2=1;
173            M1=N2;
174            M0=N1*N2;
175         }
176      }
177    
178    NFaceElements=0;    NFaceElements=0;
179    if (!periodic[0]) {    if (!periodic[0]) {
180        NDOF0=N0;        NDOF0=N0;
# Line 69  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_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
224      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
225      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, NDOF0*NDOF1*NDOF2, 0, 0 );
226      Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, NE0*NE1*NE2, NE0*NE1*NE2);
227      Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements );
228    #else
229    out=Finley_Mesh_alloc(name,3,order);    out=Finley_Mesh_alloc(name,3,order);
230    
231    if (! Finley_noError()) return NULL;    if (! Finley_noError()) return NULL;
232    
233    out->Elements=Finley_ElementFile_alloc(Hex20,out->order);    out->Elements=Finley_ElementFile_alloc(Hex20,out->order);
# Line 88  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 246  Finley_Mesh* Finley_RectangularMesh_Hex2
246    
247        
248    /*  allocate tables: */    /*  allocate tables: */
     
249    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);
250    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);
251    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
252    #endif
253    if (! Finley_noError()) {    if (! Finley_noError()) {
254        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
255        return NULL;        return NULL;
256    }    }
257      
258    #pragma omp parallel for private(i0,i1,i2,k)      /* create nodes */
259      #pragma omp parallel for private(i0,i1,i2,k)
260    for (i2=0;i2<N2;i2++) {    for (i2=0;i2<N2;i2++) {
261      for (i1=0;i1<N1;i1++) {      for (i1=0;i1<N1;i1++) {
262        for (i0=0;i0<N0;i0++) {        for (i0=0;i0<N0;i0++) {
263          k=i0+N0*i1+N0*N1*i2;          k=M0*i0+M1*i1+M2*i2;
264          out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];          out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];
265          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
266          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
267          out->Nodes->Id[k]=k;          out->Nodes->Id[k]=i0+N0*i1+N0*N1*i2;
268          out->Nodes->Tag[k]=0;          out->Nodes->Tag[k]=0;
269          out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1) +N0*N1*(i2%NDOF2);          out->Nodes->degreeOfFreedom[k]=M0*(i0%NDOF0) +M1*(i1%NDOF1) +M2*(i2%NDOF2);
270    #ifdef PASO_MPI
271            out->Nodes->Dom[k]=NODE_INTERNAL;
272    #endif
273        }        }
274      }      }
275    }    }
276    
277      
278    /* tags for the faces: */    /* tags for the faces: */
279    if (!periodic[2]) {    if (!periodic[2]) {
280       for (i1=0;i1<N1;i1++) {       for (i1=0;i1<N1;i1++) {
281         for (i0=0;i0<N0;i0++) {         for (i0=0;i0<N0;i0++) {
282           out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;           out->Nodes->Tag[M0*i0+M1*i1+M2*0]+=100;
283           out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;           out->Nodes->Tag[M0*i0+M1*i1+M2*(N2-1)]+=200;
284         }         }
285       }       }
286    }    }
287    if (!periodic[1]) {    if (!periodic[1]) {
288      for (i2=0;i2<N2;i2++) {      for (i2=0;i2<N2;i2++) {
289        for (i0=0;i0<N0;i0++) {        for (i0=0;i0<N0;i0++) {
290           out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;           out->Nodes->Tag[M0*i0+M1*0+M2*i2]+=10;
291           out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;           out->Nodes->Tag[M0*i0+M1*(N1-1)+M2*i2]+=20;
292        }        }
293      }      }
294    }    }
295    if (!periodic[0]) {    if (!periodic[0]) {
296      for (i2=0;i2<N2;i2++) {      for (i2=0;i2<N2;i2++) {
297        for (i1=0;i1<N1;i1++) {        for (i1=0;i1<N1;i1++) {
298          out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;          out->Nodes->Tag[M0*0+M1*i1+M2*i2]+=1;
299          out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;          out->Nodes->Tag[M0*(N0-1)+M1*i1+M2*i2]+=2;
300        }        }
301      }      }
302    }    }
# Line 150  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 313  Finley_Mesh* Finley_RectangularMesh_Hex2
313          out->Elements->Id[k]=k;          out->Elements->Id[k]=k;
314          out->Elements->Tag[k]=0;          out->Elements->Tag[k]=0;
315          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);
316    #ifdef PASO_MPI
317            out->Elements->Dom[k]=ELEMENT_INTERNAL;
318    #endif
319    
320          out->Elements->Nodes[INDEX2(0,k,20)]=node0;          out->Elements->Nodes[INDEX2(0,k,20)]=node0;
321          out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;          out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
# Line 200  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 366  Finley_Mesh* Finley_RectangularMesh_Hex2
366          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
367          out->FaceElements->Tag[k]=100;          out->FaceElements->Tag[k]=100;
368          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
369    #ifdef PASO_MPI
370            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
371    #endif
372        
373          if  (useElementsOnFace) {          if  (useElementsOnFace) {
374             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 247  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 416  Finley_Mesh* Finley_RectangularMesh_Hex2
416          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;
417          out->FaceElements->Tag[k]=200;          out->FaceElements->Tag[k]=200;
418          out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;          out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
419    #ifdef PASO_MPI
420            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
421    #endif
422    
423          if  (useElementsOnFace) {          if  (useElementsOnFace) {
424             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;
# Line 301  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 473  Finley_Mesh* Finley_RectangularMesh_Hex2
473           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
474           out->FaceElements->Tag[k]=1;           out->FaceElements->Tag[k]=1;
475           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
476    #ifdef PASO_MPI
477            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
478    #endif
479    
480           if  (useElementsOnFace) {           if  (useElementsOnFace) {
481              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 354  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 529  Finley_Mesh* Finley_RectangularMesh_Hex2
529           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;
530           out->FaceElements->Tag[k]=2;           out->FaceElements->Tag[k]=2;
531           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
532    #ifdef PASO_MPI
533            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
534    #endif
535    
536           if  (useElementsOnFace) {           if  (useElementsOnFace) {
537              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;
# Line 409  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 587  Finley_Mesh* Finley_RectangularMesh_Hex2
587           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
588           out->FaceElements->Tag[k]=10;           out->FaceElements->Tag[k]=10;
589           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;
590    #ifdef PASO_MPI
591            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
592    #endif
593    
594           if  (useElementsOnFace) {           if  (useElementsOnFace) {
595              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;
# Line 462  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 643  Finley_Mesh* Finley_RectangularMesh_Hex2
643           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;
644           out->FaceElements->Tag[k]=20;           out->FaceElements->Tag[k]=20;
645           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;
646    #ifdef PASO_MPI
647            out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
648    #endif
649    
650           if  (useElementsOnFace) {           if  (useElementsOnFace) {
651              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;
# Line 506  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 690  Finley_Mesh* Finley_RectangularMesh_Hex2
690    out->FaceElements->minColor=0;    out->FaceElements->minColor=0;
691    out->FaceElements->maxColor=24;    out->FaceElements->maxColor=24;
692    
693    /*  face elements done: */  #ifdef PASO_MPI
694          Finley_ElementFile_setDomainFlags( out->Elements );
695        Finley_ElementFile_setDomainFlags( out->FaceElements );
696        Finley_ElementFile_setDomainFlags( out->ContactElements );
697        Finley_ElementFile_setDomainFlags( out->Points );
698    
699        /* reorder the degrees of freedom */
700        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
701    #endif
702        
703    /*   condense the nodes: */    /*   condense the nodes: */
     
704    Finley_Mesh_resolveNodeIds(out);    Finley_Mesh_resolveNodeIds(out);
705    
706    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
# Line 518  Finley_Mesh* Finley_RectangularMesh_Hex2 Line 709  Finley_Mesh* Finley_RectangularMesh_Hex2
709    #ifdef Finley_TRACE    #ifdef Finley_TRACE
710    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
711    #endif    #endif
712      if (Finley_noError()) {
713           if (! Finley_Mesh_isPrepared(out) ) {
714              Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
715           }
716      } else {
717          Finley_Mesh_dealloc(out);
718      }
719      return out;
720    }
721    #ifdef PASO_MPI
722    Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {
723      dim_t N0,N0t,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,totalNECount,faceNECount,NDOF0,NDOF0t,NDOF1,NDOF2,NFaceElements,NUMNODES,M0,M1,M2;
724      dim_t kk,iI, NE0_local, numNodesLocal, numDOFLocal, numElementsLocal, numElementsInternal, nodesExternal[2], DOFExternal[2], numNodesExternal, DOFBoundary[2];
725        bool_t dom_left, dom_right, dom_internal;
726      index_t firstNode=0, DOFcount=0, node0, node1, node2, idCount;
727      index_t targetDomain=-1, firstNodeConstruct, j;
728      bool_t periodicLocal[2], domLeft=FALSE, domRight=FALSE, domInternal=FALSE, boundaryLeft=FALSE, boundaryRight=FALSE;
729        index_t *indexBackward=NULL, *indexForward=NULL,*facePerm=NULL, *forwardDOF=NULL, *backwardDOF=NULL;
730      Finley_Mesh* out;
731      char name[50];
732      double time0=Finley_timer();
733      Paso_MPIInfo *mpi_info = NULL;
734    
735      NE0=MAX(1,numElements[0]);
736      NE1=MAX(1,numElements[1]);
737      NE2=MAX(1,numElements[2]);
738      N0=2*NE0+1;
739      N1=2*NE1+1;
740      N2=2*NE2+1;
741    
742        index_t face0[] = {0,4,7,3,1,5,6,2,12,19,15,11,8,16,18,10,13,17,14,9};
743        index_t face1[] = {1,2,6,5,0,3,7,4,9,14,17,13,8,10,18,16,11,15,19,12};
744        index_t face2[] = {0,1,5,4,3,2,6,7,8,13,16,12,11,9,17,19,14,18,15,10};
745        index_t face3[] = {3,7,6,2,0,4,5,1,15,18,14,10,11,19,17,9,12,16,13,8};
746        index_t face4[] = {0,3,2,1,4,7,6,5,11,10,9,8,15,14,13,12,19,18,17,16};
747        index_t face5[] = {4,5,6,7,0,1,2,3,16,17,18,19,12,13,14,15,8,9,10,11};
748        
749        index_t face0R[] = {0,4,7,3,12,19,15,11};
750        index_t face1R[] = {1,2,6,5,9,14,17,13};
751        index_t face2R[] = {0,1,5,4,8,13,16,12};
752        index_t face3R[] = {3,7,6,2,15,18,14,10};
753        index_t face4R[] = {0,3,2,1,11,10,9,8};
754        index_t face5R[] = {4,5,6,7,16,17,18,19};
755        
756      /* get MPI information */
757      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
758      if (! Finley_noError())
759            return NULL;
760    
761      /* use the serial version to generate the mesh for the 1-CPU case */
762      if( mpi_info->size==1 )
763      {
764        out =  Finley_RectangularMesh_Hex20_singleCPU( numElements, Length, periodic, order, useElementsOnFace, mpi_info);
765            return out;
766      }    
767    
768      if( mpi_info->rank==0 )
769        domLeft = TRUE;
770      if( mpi_info->rank==mpi_info->size-1 )
771        domRight = TRUE;
772      if( mpi_info->rank>0 && mpi_info->rank<mpi_info->size-1 )
773        domInternal = TRUE;
774    
775      /* dimensions of the local subdomain */
776      domain_calculateDimension( mpi_info->rank, mpi_info->size, NE0, periodic[0], &numNodesLocal, &numDOFLocal, &numElementsLocal, &numElementsInternal, &firstNode, nodesExternal, DOFExternal, &numNodesExternal, periodicLocal, DOFBoundary );  
777    
778      NFaceElements=0;
779      if (!periodic[0]) {
780          NDOF0=N0;
781          NFaceElements+=(domRight+domLeft)*NE1*NE2;
782      } else {
783          NDOF0=N0-1;
784      }
785      if (!periodic[1]) {
786          NDOF1=N1;
787          NFaceElements+=2*numElementsLocal*NE2;
788      } else {
789          NDOF1=N1-1;
790      }
791      if (!periodic[2]) {
792          NDOF2=N2;
793          NFaceElements+=2*numElementsLocal*NE1;
794      } else {
795          NDOF2=N2-1;
796      }
797      
798        boundaryLeft = !domLeft || periodicLocal[0];
799        boundaryRight = !domRight || periodicLocal[1];
800        N0t = numNodesLocal + boundaryRight + boundaryLeft*2;
801        NDOF0t = numDOFLocal + boundaryRight + boundaryLeft*2;
802        firstNodeConstruct = firstNode - 2*boundaryLeft;
803        firstNodeConstruct = firstNodeConstruct<0 ? N0+firstNodeConstruct-1 : firstNodeConstruct;
804    
805      /*  allocate mesh: */
806      
807      sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
808      out=Finley_Mesh_alloc(name,3,order,mpi_info);
809    
810      if (! Finley_noError()) return NULL;
811    
812      out->Elements=Finley_ElementFile_alloc(Hex20,out->order,mpi_info);
813      if (useElementsOnFace) {
814         out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order,mpi_info);
815         out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order,mpi_info);
816      } else {
817         out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order,mpi_info);
818         out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order,mpi_info);
819      }
820      out->Points=Finley_ElementFile_alloc(Point1,out->order,mpi_info);
821      if (! Finley_noError()) {
822          Finley_Mesh_dealloc(out);
823          return NULL;
824      }
825    
826      /*  allocate tables: */
827      Finley_NodeFile_allocTable(out->Nodes,N0t*N1*N2);
828      Finley_ElementFile_allocTable(out->Elements,numElementsLocal*NE1*NE2);
829      Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
830        
831      Finley_NodeDistribution_allocTable( out->Nodes->degreeOfFreedomDistribution, numDOFLocal*NDOF1*NDOF2, NDOF1*NDOF2*3, 0 );
832      Finley_ElementDistribution_allocTable( out->Elements->elementDistribution, numElementsLocal*NE1*NE2, NE1*NE2*(numElementsLocal-boundaryRight*(!periodic[1])) );
833      Finley_ElementDistribution_allocTable( out->FaceElements->elementDistribution, NFaceElements, NFaceElements-2*boundaryRight*(NE2*(!periodic[1])+NE1*(!periodic[2])) );
834    if (! Finley_noError()) {    if (! Finley_noError()) {
835        Finley_Mesh_dealloc(out);        Finley_Mesh_dealloc(out);
836        return NULL;        return NULL;
837    }    }
   return out;  
 }  
838    
839  /*      k=0;
840  * Revision 1.3  2005/09/01 03:31:35  jgs    #pragma omp parallel for private(i0,i1,i2,k)
841  * Merge of development branch dev-02 back to main trunk on 2005-09-01    for (i2=0;i2<N2;i2++) {
842  *      for (i1=0;i1<N1;i1++) {
843  * Revision 1.2.2.2  2005/09/07 06:26:19  gross        for (i0=0;i0<N0t;i0++,k++) {        
844  * the solver from finley are put into the standalone package paso now          out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE((i0+firstNodeConstruct) % N0)/DBLE(N0-1)*Length[0];
845  *          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];
846  * Revision 1.2.2.1  2005/08/24 02:02:18  gross          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];
847  * timing output switched off. solver output can be swiched through getSolution(verbose=True) now.          out->Nodes->Id[k]=k;
848  *          out->Nodes->Tag[k]=0;
849  * Revision 1.2  2005/07/08 04:07:52  jgs          out->Nodes->degreeOfFreedom[k]=i0 + (i1%NDOF1)*N0t + (i2%NDOF2)*N0t*N1;
850  * Merge of development branch back to main trunk on 2005-07-08                  out->Nodes->Dom[k]=NODE_INTERNAL;
851  *        }
852  * Revision 1.1.1.1.2.1  2005/06/29 02:34:51  gross      }
853  * some changes towards 64 integers in finley    }
854  *  
855  * Revision 1.1.1.1  2004/10/26 06:53:57  jgs      /* mark the nodes that reference external and boundary DOF as such */
856  * initial import of project esys2      if( boundaryLeft ){
857  *          for( i1=0; i1<N1; i1++ )
858  * Revision 1.1.1.1  2004/06/24 04:00:40  johng              for( i2=0; i2<N2; i2++ ) {
859  * Initial version of eys using boost-python.                  out->Nodes->Dom[N1*N0t*i2+N0t*i1] = NODE_EXTERNAL;
860  *                  out->Nodes->Dom[N1*N0t*i2+N0t*i1+1] = NODE_EXTERNAL;
861  *                  out->Nodes->Dom[N1*N0t*i2+N0t*i1+2] = NODE_BOUNDARY;
862  */              }
863        }
864        if( boundaryRight ){
865            for( i1=0; i1<N1; i1++ )
866                for( i2=0; i2<N2; i2++ ) {
867                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-1] = NODE_EXTERNAL;
868                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-2] = NODE_BOUNDARY;
869                    out->Nodes->Dom[N1*N0t*i2+N0t*(i1+1)-3] = NODE_BOUNDARY;
870                }
871        }
872        if( periodicLocal[0] ){
873            for( i1=0; i1<N1; i1++ )
874                for( i2=0; i2<N2; i2++ ) {
875                    out->Nodes->degreeOfFreedom[N1*N0t*i2+i1*N0t+3] = out->Nodes->degreeOfFreedom[i1*N0t+2];
876                    out->Nodes->Dom[N1*N0t*i2+N0t*i1+3] = NODE_BOUNDARY;
877                }
878        }
879            
880      /* tag Nodes that are referenced by face elements */
881      if (!periodic[2]) {    
882        for (i1=0;i1<N1;i1++) {
883          for (i0=0;i0<N0t;i0++) {  
884             out->Nodes->Tag[i0 + N0t*i1]+=100;
885             out->Nodes->Tag[i0 + N0t*i1 + N0t*N1*(N2-1)]+=200;
886           }
887         }
888      }
889      if (!periodic[1]) {
890        for (i2=0;i2<N2;i2++) {
891          for (i0=0;i0<N0t;i0++) {
892             out->Nodes->Tag[i0 + i2*N1*N0t]+=10;
893             out->Nodes->Tag[i0 + (i2+1)*N1*N0t-N0t]+=20;
894          }
895        }
896      }
897      if (!periodic[0] && !domInternal ) {
898        for (i2=0;i2<N2;i2++) {
899          for (i1=0;i1<N1;i1++) {
900            if( domLeft )
901              out->Nodes->Tag[i1*N0t + i2*N0t*N1]+=1;
902            if( domRight )
903              out->Nodes->Tag[(i1+1)*N0t-1 + i2*N0t*N1]+=2;
904          }
905        }
906      }
907    
908        /* form the boudary communication information */
909        forwardDOF  = MEMALLOC(NDOF1*NDOF2*2,index_t);
910        backwardDOF = MEMALLOC(NDOF1*NDOF2*2,index_t);
911        if( !(mpi_info->size==2 && periodicLocal[0])){
912            if( boundaryLeft  ) {
913                targetDomain = mpi_info->rank-1 < 0 ? mpi_info->size-1 : mpi_info->rank-1;
914                for( iI=0, i2=0; i2<NDOF2; i2++ ){
915                    for( i1=0; i1<NDOF1; i1++, iI+=2 ){
916                        forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
917                        backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
918                        backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
919                    }
920                }
921                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
922                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
923            }
924            if( boundaryRight ) {
925                targetDomain = mpi_info->rank+1 > mpi_info->size-1 ? 0 : mpi_info->rank+1;
926                for( iI=0,i2=0; i2<NDOF2; i2++ ){
927                    for( i1=0; i1<NDOF1; i1++, iI+=2 ){
928                        forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
929                        forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
930                        backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
931                    }
932                }
933                Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
934                Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
935            }
936        } else{
937            /* periodic boundary conditions with 2 domains, need to change the order in which domain 0 passes boundary data */
938            targetDomain = 1;
939            
940            for( iI=0,i2=0; i2<NDOF2; i2++ ){
941                for( i1=0; i1<NDOF1; i1++, iI+=2 ){
942                    forwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-3];
943                    forwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-2];
944                    backwardDOF[i1+i2*NDOF1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+(i1+1)*N0t-1];
945                }
946            }
947            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, forwardDOF );
948            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, backwardDOF );
949    
950            for( iI=0, i2=0; i2<NDOF2; i2++ ){
951                for( i1=0; i1<NDOF1; i1++, iI+=2 ){
952                    forwardDOF[i1+i2*NDOF1]  = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+2];
953                    backwardDOF[iI] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t];
954                    backwardDOF[iI+1] = out->Nodes->degreeOfFreedom[i2*N0t*N1+i1*N0t+1];
955                }
956            }
957            Finley_NodeDistribution_addForward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2, forwardDOF );
958            Finley_NodeDistribution_addBackward( out->Nodes->degreeOfFreedomDistribution, targetDomain, NDOF1*NDOF2*2, backwardDOF );
959        }
960        MEMFREE( forwardDOF );
961        MEMFREE( backwardDOF );
962      /*   set the elements: */
963      
964        k=0;
965      #pragma omp parallel for private(i0,i1,i2,k,node0)
966      for (i2=0;i2<NE2;i2++) {
967        for (i1=0;i1<NE1;i1++) {
968          for (i0=0;i0<numElementsLocal;i0++,k++) {
969                    node0 = (periodicLocal[0] && !i0) ? 2*(i1*N0t + i2*N1*N0t) : 2*(i1*N0t + i2*N1*N0t + i0) + periodicLocal[0];
970    
971                    out->Elements->Id[k]=((firstNodeConstruct/2+i0)%NE0)*NE1*NE2 + NE1*i2 + i1;
972            out->Elements->Tag[k]=0;
973            out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);
974            out->Elements->Dom[k]=ELEMENT_INTERNAL;
975    
976            out->Elements->Nodes[INDEX2(0,k,20)]=node0;
977            out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;
978            out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0t+2;
979            out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0t;
980            out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0t*N1;
981            out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0t*N1+2;
982            out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0t*N1+2*N0t+2;
983            out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0t*N1+2*N0t;
984            out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;
985            out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0t+2;
986            out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0t+1;
987            out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0t;
988            out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0t*N1;
989            out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0t*N1+2;
990            out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0t*N1+2*N0t+2;
991            out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0t*N1+2*N0t;
992            out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0t*N1+1;
993            out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0t*N1+N0t+2;
994            out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0t*N1+2*N0t+1;
995            out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0t*N1+N0t;
996          }
997        }
998      }
999      out->Elements->minColor=0;
1000      out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);
1001    
1002        if( boundaryLeft )
1003            for( i2=0; i2<NE2; i2++ )
1004                for( i1=0; i1<NE1; i1++ )
1005                    out->Elements->Dom[i2*NE1*numElementsLocal+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1006        if( boundaryRight )
1007            for( i2=0; i2<NE2; i2++ )
1008                for( i1=0; i1<NE1; i1++ )
1009                    out->Elements->Dom[i2*NE1*numElementsLocal+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1010    
1011        out->Elements->numElements = numElementsLocal*NE1*NE2;
1012        Finley_ElementFile_setDomainFlags( out->Elements );
1013      
1014      /*   face elements: */
1015      
1016      if  (useElementsOnFace) {
1017         NUMNODES=20;
1018      } else {
1019         NUMNODES=8;
1020      }
1021      totalNECount=out->Elements->numElements;
1022      faceNECount=0;
1023        idCount = totalNECount;
1024      
1025      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
1026        numElementsInternal = numElementsLocal-nodesExternal[0]-nodesExternal[1];
1027      if (!periodic[0] && !domInternal) {
1028         /* **  elements on boundary 001 (x1=0): */
1029         if( domLeft ){
1030                 #pragma omp parallel for private(i1,i2,k)
1031                 for (i2=0;i2<NE2;i2++) {
1032                     for (i1=0;i1<NE1;i1++) {
1033                         k=i1+NE1*i2+faceNECount;
1034                         kk=i1*numElementsLocal + i2*numElementsLocal*NE1;
1035                         facePerm =!useElementsOnFace ? face0R :  face0;
1036            
1037                         out->FaceElements->Id[k]=idCount++;
1038                         out->FaceElements->Tag[k]=1;
1039                         out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1040                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;
1041    
1042                         for( j=0; j<NUMNODES; j++ )
1043                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1044                     }
1045                 }
1046                 totalNECount+=NE1*NE2;
1047                 faceNECount+=NE1*NE2;
1048         }
1049         /* **  elements on boundary 002 (x1=1): */
1050             if( domRight ) {
1051                 #pragma omp parallel for private(i1,i2,k)
1052                 for (i2=0;i2<NE2;i2++) {
1053                     for (i1=0;i1<NE1;i1++) {
1054                         k=i1+NE1*i2+faceNECount;
1055                         kk=(i1+1)*numElementsLocal + i2*numElementsLocal*NE1 - 1;
1056                         facePerm =!useElementsOnFace ? face1R :  face1;
1057            
1058                         out->FaceElements->Id[k]=idCount++;
1059                         out->FaceElements->Tag[k]=2;
1060                 out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1061                         out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;
1062    
1063                         for( j=0; j<NUMNODES; j++ )
1064                            out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1065                     }
1066                 }
1067           totalNECount+=NE1*NE2;
1068           faceNECount+=NE1*NE2;
1069         }
1070      }
1071      if (!periodic[1]) {
1072         /* **  elements on boundary 010 (x2=0): */
1073      
1074         #pragma omp parallel for private(i0,i2,k)
1075         for (i2=0;i2<NE2;i2++) {
1076           for (i0=0;i0<numElementsLocal;i0++) {
1077             k=i0+numElementsLocal*i2+faceNECount;
1078                     kk=i0+numElementsLocal*NE1*i2;
1079                     facePerm =!useElementsOnFace ? face2R :  face2;
1080      
1081             out->FaceElements->Id[k]=idCount++;
1082             out->FaceElements->Tag[k]=10;
1083             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1084             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;
1085    
1086                     for( j=0; j<NUMNODES; j++ )
1087                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1088           }
1089         }
1090           if( boundaryLeft ){
1091                for( i2=0; i2<NE2; i2++ )
1092                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1093                if( periodicLocal[0] )
1094                    for( i2=0; i2<NE2; i2++ )
1095                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1096             }
1097           if( boundaryRight )
1098                for( i2=0; i2<NE2; i2++ )
1099                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1100         totalNECount+=numElementsLocal*NE2;
1101         faceNECount+=numElementsLocal*NE2;
1102      
1103         /* **  elements on boundary 020 (x2=1): */
1104      
1105         #pragma omp parallel for private(i0,i2,k)
1106         for (i2=0;i2<NE2;i2++) {
1107           for (i0=0;i0<numElementsLocal;i0++) {
1108             k=i0+numElementsLocal*i2+faceNECount;
1109                     kk=i0+numElementsLocal*NE1*(i2+1)-numElementsLocal;
1110                     facePerm =!useElementsOnFace ? face3R :  face3;
1111      
1112             out->FaceElements->Tag[k]=20;
1113             out->FaceElements->Id[k]=idCount++;
1114             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1115             out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;
1116    
1117                     for( j=0; j<NUMNODES; j++ )
1118                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1119           }
1120         }
1121           if( boundaryLeft ){
1122                for( i2=0; i2<NE2; i2++ )
1123                    out->FaceElements->Dom[faceNECount+i2*numElementsLocal]=ELEMENT_BOUNDARY;
1124                if( periodicLocal[0] )
1125                    for( i2=0; i2<NE2; i2++ )
1126                        out->FaceElements->Dom[faceNECount+i2*numElementsLocal+1]=ELEMENT_BOUNDARY;
1127             }
1128           if( boundaryRight )
1129                for( i2=0; i2<NE2; i2++ )
1130                    out->FaceElements->Dom[faceNECount+(i2+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1131         totalNECount+=numElementsLocal*NE2;
1132         faceNECount+=numElementsLocal*NE2;
1133      }
1134      if (!periodic[2]) {
1135         /*  elements on boundary 100 (x3=0): */
1136      
1137         #pragma omp parallel for private(i0,i1,k)
1138         for (i1=0;i1<NE1;i1++) {
1139           for (i0=0; i0<numElementsLocal; i0++) {
1140             k=i0+numElementsLocal*i1+faceNECount;
1141                     kk=i0 + i1*numElementsLocal;
1142                     facePerm =!useElementsOnFace ? face4R :  face4;
1143      
1144             out->FaceElements->Id[k]=idCount++;
1145             out->FaceElements->Tag[k]=100;
1146             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1147             out->FaceElements->Color[k]=(i0%2)+2*(i1%2);
1148    
1149                     for( j=0; j<NUMNODES; j++ )
1150                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1151           }
1152         }
1153           if( boundaryLeft ){
1154                for( i1=0; i1<NE1; i1++ )
1155                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1156                if( periodicLocal[0] )
1157                    for( i1=0; i1<NE1; i1++ )
1158                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1159             }
1160           if( boundaryRight )
1161                for( i1=0; i1<NE1; i1++ )
1162                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1163         totalNECount+=NE1*numElementsLocal;
1164         faceNECount+=NE1*numElementsLocal;
1165            
1166         /* **  elements on boundary 200 (x3=1) */
1167      
1168         #pragma omp parallel for private(i0,i1,k)
1169         for (i1=0;i1<NE1;i1++) {
1170           for (i0=0;i0<numElementsLocal;i0++) {
1171             k=i0+numElementsLocal*i1+faceNECount;
1172                     kk=i0+i1*numElementsLocal+numElementsLocal*NE1*(NE2-1);
1173                     facePerm = !useElementsOnFace ? face5R : face5;
1174    
1175             out->FaceElements->Id[k]=idCount++;
1176             out->FaceElements->Tag[k]=200;
1177             out->FaceElements->Dom[k]=ELEMENT_INTERNAL;
1178             out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;
1179    
1180                     for( j=0; j<NUMNODES; j++ )
1181                        out->FaceElements->Nodes[INDEX2(j,k,NUMNODES)]=out->Elements->Nodes[INDEX2(facePerm[j],kk,20)];
1182           }
1183         }
1184           if( boundaryLeft ){
1185                for( i1=0; i1<NE1; i1++ )
1186                    out->FaceElements->Dom[faceNECount+i1*numElementsLocal]=ELEMENT_BOUNDARY;
1187                if( periodicLocal[0] )
1188                    for( i1=0; i1<NE1; i1++ )
1189                        out->FaceElements->Dom[faceNECount+i1*numElementsLocal+1]=ELEMENT_BOUNDARY;
1190             }
1191           if( boundaryRight )
1192                for( i1=0; i1<NE1; i1++ )
1193                    out->FaceElements->Dom[faceNECount+(i1+1)*numElementsLocal-1]=ELEMENT_BOUNDARY;
1194         totalNECount+=NE1*numElementsLocal;
1195         faceNECount+=NE1*numElementsLocal;
1196      }
1197        out->FaceElements->elementDistribution->numInternal = faceNECount;
1198        
1199      out->FaceElements->minColor=0;
1200      out->FaceElements->maxColor=23;
1201        out->FaceElements->numElements=faceNECount;
1202        Finley_ElementFile_setDomainFlags( out->FaceElements );
1203    
1204      /* setup distribution info for other elements */
1205        Finley_ElementFile_setDomainFlags( out->ContactElements );
1206        Finley_ElementFile_setDomainFlags( out->Points );
1207    
1208        /* reorder the degrees of freedom */
1209        Finley_Mesh_resolveDegreeOfFreedomOrder( out, TRUE );
1210    
1211      /*   condense the nodes: */
1212      Finley_Mesh_resolveNodeIds(out);
1213      if( !Finley_MPI_noError(mpi_info) )
1214      {
1215        Paso_MPIInfo_dealloc( mpi_info );
1216        Finley_Mesh_dealloc(out);
1217        return NULL;
1218      }
1219    
1220      /* prepare mesh for further calculatuions:*/
1221      Finley_Mesh_prepare(out);
1222      if( !Finley_MPI_noError(mpi_info) )
1223      {
1224        Paso_MPIInfo_dealloc( mpi_info );
1225        Finley_Mesh_dealloc(out);
1226        return NULL;
1227      }
1228    
1229      /* free up memory */
1230      Paso_MPIInfo_dealloc( mpi_info );
1231    
1232        //print_mesh_statistics( out, FALSE );
1233    
1234      #ifdef Finley_TRACE
1235      printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
1236      #endif
1237      if (Finley_noError()) {
1238           if (!Finley_Mesh_isPrepared(out) ) {
1239              Finley_setError(SYSTEM_ERROR,"Mesh is not prepared for calculation. Contact the programmers.");
1240           }
1241      }
1242      return out;
1243    }
1244    #endif
1245    

Legend:
Removed from v.150  
changed lines
  Added in v.964

  ViewVC Help
Powered by ViewVC 1.1.26