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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26