/[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

trunk/esys2/finley/src/finleyC/Mesh_hex8.c revision 82 by jgs, Tue Oct 26 06:53:54 2004 UTC trunk/finley/src/Mesh_hex8.c revision 2881 by jfenwick, Thu Jan 28 02:03:15 2010 UTC
# Line 1  Line 1 
1    
2    /*******************************************************
3    *
4    * Copyright (c) 2003-2010 by University of Queensland
5    * Earth Systems Science Computational Center (ESSCC)
6    * http://www.uq.edu.au/esscc
7    *
8    * Primary Business: Queensland, Australia
9    * Licensed under the Open Software License version 3.0
10    * http://www.opensource.org/licenses/osl-3.0.php
11    *
12    *******************************************************/
13    
14    
15  /**************************************************************/  /**************************************************************/
16    
17  /*   Finley: generates rectangular meshes */  /*   Finley: generates rectangular meshes */
# Line 8  Line 22 
22    
23  /**************************************************************/  /**************************************************************/
24    
 /*   Copyrights by ACcESS Australia 2003/04 */  
 /*   Author: gross@access.edu.au */  
 /*   Version: $Id$ */  
   
 /**************************************************************/  
   
 #include "Common.h"  
 #include "Finley.h"  
 #include "Mesh.h"  
25  #include "RectangularMesh.h"  #include "RectangularMesh.h"
26    
27  /**************************************************************/  Finley_Mesh* Finley_RectangularMesh_Hex8(dim_t* numElements,
28                                              double* Length,
29  Finley_Mesh* Finley_RectangularMesh_Hex8(int* numElements,double* Length,int* periodic, int order,int useElementsOnFace) {                                            bool_t* periodic,
30    int N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,node0,totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements,NUMNODES;                                            index_t order,
31                                              index_t reduced_order,
32                                              bool_t useElementsOnFace,
33                                              bool_t useFullElementOrder,
34                                              bool_t optimize)
35    {
36      #define N_PER_E 1
37      #define DIM 3
38      dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0=0, Nstride1=0,Nstride2=0, local_NE0, local_NE1, local_NE2, local_N0=0, local_N1=0, local_N2=0;
39      dim_t totalNECount,faceNECount,NDOF0=0,NDOF1=0,NDOF2=0,NFaceElements=0, NN;
40      index_t node0, myRank, e_offset2, e_offset1, e_offset0=0, offset1=0, offset2=0, offset0=0, global_i0, global_i1, global_i2;
41      Finley_ReferenceElementSet *refPoints=NULL, *refContactElements=NULL, *refFaceElements=NULL, *refElements=NULL;
42    Finley_Mesh* out;    Finley_Mesh* out;
43      Paso_MPIInfo *mpi_info = NULL;
44    char name[50];    char name[50];
45      #ifdef Finley_TRACE
46    double time0=Finley_timer();    double time0=Finley_timer();
47      #endif
48    
49      /* get MPI information */
50      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
51      if (! Finley_noError()) {
52            return NULL;
53      }
54      myRank=mpi_info->rank;
55    
56      /* set up the global dimensions of the mesh */
57    
58    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
59    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
60    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
61    N0=NE0+1;    N0=N_PER_E*NE0+1;
62    N1=NE1+1;    N1=N_PER_E*NE1+1;
63    N2=NE2+1;    N2=N_PER_E*NE2+1;
   
   NFaceElements=0;  
   if (!periodic[0]) {  
       NDOF0=N0;  
       NFaceElements+=2*NE1*NE2;  
   } else {  
       NDOF0=N0-1;  
   }  
   if (!periodic[1]) {  
       NDOF1=N1;  
       NFaceElements+=2*NE0*NE2;  
   } else {  
       NDOF1=N1-1;  
   }  
   if (!periodic[2]) {  
       NDOF2=N2;  
       NFaceElements+=2*NE0*NE1;  
   } else {  
       NDOF2=N2-1;  
   }  
     
   /*  allocate mesh: */  
     
   sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);  
   out=Finley_Mesh_alloc(name,3,order);  
   if (Finley_ErrorCode!=NO_ERROR) return NULL;  
64    
65    out->Elements=Finley_ElementFile_alloc(Hex8,out->order);    /*  allocate mesh: */  
66    if (useElementsOnFace) {    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);
67       out->FaceElements=Finley_ElementFile_alloc(Hex8Face,out->order);    out=Finley_Mesh_alloc(name,DIM, mpi_info);
68       out->ContactElements=Finley_ElementFile_alloc(Hex8Face_Contact,out->order);    if (! Finley_noError()) {
69    } else {        Paso_MPIInfo_free( mpi_info );
      out->FaceElements=Finley_ElementFile_alloc(Rec4,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Rec4_Contact,out->order);  
   }  
   out->Points=Finley_ElementFile_alloc(Point1,out->order);  
   if (Finley_ErrorCode!=NO_ERROR) {  
       Finley_Mesh_dealloc(out);  
70        return NULL;        return NULL;
71    }    }
72      refElements= Finley_ReferenceElementSet_alloc(Hex8,order,reduced_order);
73        if (useElementsOnFace) {
74    /*  allocate tables: */          refFaceElements=Finley_ReferenceElementSet_alloc(Hex8Face, order, reduced_order);
75              refContactElements=Finley_ReferenceElementSet_alloc(Hex8Face_Contact, order, reduced_order);
76    Finley_NodeFile_allocTable(out->Nodes,N0*N1*N2);    } else {
77    Finley_ElementFile_allocTable(out->Elements,NE0*NE1*NE2);          refFaceElements=Finley_ReferenceElementSet_alloc(Rec4, order, reduced_order);
78    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);          refContactElements=Finley_ReferenceElementSet_alloc(Rec4_Contact, order, reduced_order);
   if (Finley_ErrorCode!=NO_ERROR) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
79    }    }
80      refPoints=Finley_ReferenceElementSet_alloc(Point1, order, reduced_order);
81        
82    /*  set nodes: */  
83      if ( Finley_noError()) {
84        
85    #pragma omp parallel for private(i0,i1,i2,k)        Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(refPoints, mpi_info));
86    for (i2=0;i2<N2;i2++) {        Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(refContactElements, mpi_info));
87      for (i1=0;i1<N1;i1++) {        Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(refFaceElements, mpi_info));
88        for (i0=0;i0<N0;i0++) {        Finley_Mesh_setElements(out,Finley_ElementFile_alloc(refElements, mpi_info));
89          k=i0+N0*i1+N0*N1*i2;  
90          out->Nodes->Coordinates[INDEX2(0,k,3)]=DBLE(i0)/DBLE(N0-1)*Length[0];        /* work out the largest dimension */
91          out->Nodes->Coordinates[INDEX2(1,k,3)]=DBLE(i1)/DBLE(N1-1)*Length[1];        if (N2==MAX3(N0,N1,N2)) {
92          out->Nodes->Coordinates[INDEX2(2,k,3)]=DBLE(i2)/DBLE(N2-1)*Length[2];            Nstride0=1;
93          out->Nodes->Id[k]=k;            Nstride1=N0;
94          out->Nodes->Tag[k]=0;            Nstride2=N0*N1;
95          out->Nodes->degreeOfFreedom[k]=(i0%NDOF0) +N0*(i1%NDOF1) +N0*N1*(i2%NDOF2);            local_NE0=NE0;
96        }            e_offset0=0;
97      }            local_NE1=NE1;
98    }            e_offset1=0;
99    /* tags for the faces: */            Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
100    if (!periodic[2]) {        } else if (N1==MAX3(N0,N1,N2)) {
101      for (i1=0;i1<N1;i1++) {            Nstride0=N2;
102        for (i0=0;i0<N0;i0++) {            Nstride1=N0*N2;
103          out->Nodes->Tag[i0+N0*i1+N0*N1*0]+=100;            Nstride2=1;
104          out->Nodes->Tag[i0+N0*i1+N0*N1*(N2-1)]+=200;            local_NE0=NE0;
105        }          e_offset0=0;
106      }          Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
107    }          local_NE2=NE2;
108    if (!periodic[1]) {          e_offset2=0;
109      for (i2=0;i2<N2;i2++) {        } else {
110        for (i0=0;i0<N0;i0++) {            Nstride0=N1*N2;
111           out->Nodes->Tag[i0+N0*0+N0*N1*i2]+=10;            Nstride1=1;
112           out->Nodes->Tag[i0+N0*(N1-1)+N0*N1*i2]+=20;            Nstride2=N1;
113        }            Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
114      }            local_NE1=NE1;
115    }            e_offset1=0;
116    if (!periodic[0]) {            local_NE2=NE2;
117      for (i2=0;i2<N2;i2++) {            e_offset2=0;
118        for (i1=0;i1<N1;i1++) {        }
119          out->Nodes->Tag[0+N0*i1+N0*N1*i2]+=1;        offset0=e_offset0*N_PER_E;
120          out->Nodes->Tag[(N0-1)+N0*i1+N0*N1*i2]+=2;        offset1=e_offset1*N_PER_E;
121        }        offset2=e_offset2*N_PER_E;
122      }        local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
123          local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
124          local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
125    
126          /* get the number of surface elements */
127    
128          NFaceElements=0;
129          if (!periodic[2] && (local_NE2>0)) {
130              NDOF2=N2;
131              if (offset2==0) NFaceElements+=local_NE1*local_NE0;
132              if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
133          } else {
134              NDOF2=N2-1;
135          }
136    
137          if (!periodic[0] && (local_NE0>0)) {
138              NDOF0=N0;
139              if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
140              if (local_NE0+e_offset0 == NE0) NFaceElements+=local_NE1*local_NE2;
141          } else {
142              NDOF0=N0-1;
143          }
144             if (!periodic[1] && (local_NE1>0)) {
145                 NDOF1=N1;
146                 if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
147                 if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
148             } else {
149                 NDOF1=N1-1;
150             }
151    }    }
152        
   /*   set the elements: */  
153        
154    #pragma omp parallel for private(i0,i1,i2,k,node0)    /*  allocate tables: */
155    for (i2=0;i2<NE2;i2++) {    if (Finley_noError()) {
     for (i1=0;i1<NE1;i1++) {  
       for (i0=0;i0<NE0;i0++) {  
         k=i0+NE0*i1+NE0*NE1*i2;  
         node0=i0+i1*N0+N0*N1*i2;  
   
         out->Elements->Id[k]=k;  
         out->Elements->Tag[k]=0;  
         out->Elements->Color[k]=COLOR_MOD(i0)+3*COLOR_MOD(i1)+9*COLOR_MOD(i2);;  
   
         out->Elements->Nodes[INDEX2(0,k,8)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,8)]=node0+1;  
         out->Elements->Nodes[INDEX2(2,k,8)]=node0+N0+1;  
         out->Elements->Nodes[INDEX2(3,k,8)]=node0+N0;  
         out->Elements->Nodes[INDEX2(4,k,8)]=node0+N0*N1;  
         out->Elements->Nodes[INDEX2(5,k,8)]=node0+N0*N1+1;  
         out->Elements->Nodes[INDEX2(6,k,8)]=node0+N0*N1+N0+1;  
         out->Elements->Nodes[INDEX2(7,k,8)]=node0+N0*N1+N0;  
156    
157        }        Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
158      }        Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
159          Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
160    }    }
   out->Elements->numColors=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0)+1;  
161        
162    /*   face elements: */    if (Finley_noError()) {
163           /* create nodes */
   if  (useElementsOnFace) {  
      NUMNODES=8;  
   } else {  
      NUMNODES=4;  
   }  
   totalNECount=NE0*NE1*NE2;  
   faceNECount=0;  
     
   /*   these are the quadrilateral elements on boundary 1 (x3=0): */  
   if (!periodic[2]) {  
      /* **  elements on boundary 100 (x3=0): */  
     
      #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+faceNECount;  
          node0=i0+i1*N0;  
     
          out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;  
          out->FaceElements->Tag[k]=100;  
          out->FaceElements->Color[k]=(i0%2)+2*(i1%2);  
   
          if  (useElementsOnFace) {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+1;  
          }  
        }  
      }  
      totalNECount+=NE1*NE0;  
      faceNECount+=NE1*NE0;  
     
      /* **  elements on boundary 200 (x3=1) */  
     
      #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+faceNECount;  
          node0=i0+i1*N0+N0*N1*(NE2-1);  
164        
165           out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;       #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
166           out->FaceElements->Tag[k]=200;       for (i2=0;i2<local_N2;i2++) {
167           out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;         for (i1=0;i1<local_N1;i1++) {
168             for (i0=0;i0<local_N0;i0++) {
169           if  (useElementsOnFace) {             k=i0+local_N0*i1+local_N0*local_N1*i2;
170              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;             global_i0=i0+offset0;
171              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;             global_i1=i1+offset1;
172              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;             global_i2=i2+offset2;
173              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
174              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;             out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
175              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+1;             out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
176              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+1;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
177              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;             out->Nodes->Tag[k]=0;
178           } else {             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
179              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+ N0 * N1;                                                 +Nstride1*(global_i1%NDOF1)
180              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+ N0 * N1+1;                                                 +Nstride2*(global_i2%NDOF2);
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+ N0 * N1+N0+1;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+ N0 * N1+N0;  
181           }           }
   
     
182         }         }
183       }       }
184       totalNECount+=NE1*NE0;       /*   set the elements: */
185       faceNECount+=NE1*NE0;       NN=out->Elements->numNodes;
186    }       #pragma omp parallel for private(i0,i1,i2,k,node0)
187    if (!periodic[0]) {       for (i2=0;i2<local_NE2;i2++) {
188       /* **  elements on boundary 001 (x1=0): */         for (i1=0;i1<local_NE1;i1++) {
189               for (i0=0;i0<local_NE0;i0++) {
190       #pragma omp parallel for private(i1,i2,k,node0)            
191       for (i2=0;i2<NE2;i2++) {             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
192         for (i1=0;i1<NE1;i1++) {             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
          k=i1+NE1*i2+faceNECount;  
          node0=i1*N0+N0*N1*i2;  
193        
194           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
195           out->FaceElements->Tag[k]=1;             out->Elements->Tag[k]=0;
196           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;             out->Elements->Owner[k]=myRank;
197    
198           if  (useElementsOnFace) {             out->Elements->Nodes[INDEX2(0,k,NN)] =node0                             ;
199              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                    Nstride0;
200              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;             out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
201              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;             out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1;
202              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                   ;
203              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;             out->Elements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2         +Nstride0;
204              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+1;             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
205              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0+1;             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2+Nstride1           ;
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0+1;  
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0;  
206           }           }
207         }         }
208       }       }
209       totalNECount+=NE1*NE2;       /* face elements */
210       faceNECount+=NE1*NE2;       NN=out->FaceElements->numNodes;
211         totalNECount=NE0*NE1*NE2;
212         faceNECount=0;
213         /*   these are the quadrilateral elements on boundary 1 (x3=0): */
214         if (!periodic[2]  && (local_NE2>0)) {
215           /* **  elements on boundary 100 (x3=0): */
216           if (e_offset2==0) {
217              #pragma omp parallel for private(i0,i1,k,node0)
218              for (i1=0;i1<local_NE1;i1++) {
219                for (i0=0;i0<local_NE0;i0++) {
220              
221                  k=i0+local_NE0*i1+faceNECount;
222                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
223            
224       /* **  elements on boundary 002 (x1=1): */                out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
225                    out->FaceElements->Tag[k]=100;
226       #pragma omp parallel for private(i1,i2,k,node0)                out->FaceElements->Owner[k]=myRank;
227       for (i2=0;i2<NE2;i2++) {          
228         for (i1=0;i1<NE1;i1++) {                if  (useElementsOnFace) {
229           k=i1+NE1*i2+faceNECount;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
230           node0=(NE0-1)+i1*N0+N0*N1*i2 ;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0         +Nstride1           ;
231                       out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0         +Nstride1+Nstride0;
232           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+          Nstride0           ;
233           out->FaceElements->Tag[k]=2;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride2                      ;
234           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride1           ;
235                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
236           if  (useElementsOnFace) {                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride2          +Nstride0;
237              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;                } else {
238              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
239              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           Nstride1           ;
240              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           Nstride1+Nstride0;
241              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                     Nstride0;
242              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0;                }
243              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+N0;              }
244              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1;            }
245           } else {            faceNECount+=local_NE1*local_NE0;
246              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+1;         }
247              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0+1;         totalNECount+=NE1*NE0;
248              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;         /* **  elements on boundary 200 (x3=1) */
249              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0*N1+1;         if (local_NE2+e_offset2 == NE2) {
250           }            #pragma omp parallel for private(i0,i1,k,node0)
251              for (i1=0;i1<local_NE1;i1++) {
252                for (i0=0;i0<local_NE0;i0++) {
253          
254                  k=i0+local_NE0*i1+faceNECount;
255                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
256            
257                  out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
258                  out->FaceElements->Tag[k]=200;
259                  out->FaceElements->Owner[k]=myRank;
260                  if  (useElementsOnFace) {
261                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
262                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2+           Nstride0;
263                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
264                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
265          
266                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
267                     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride0                      ;
268                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           Nstride1+Nstride0;
269                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           Nstride1;
270                  } else {
271                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+Nstride2                      ;
272                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2           +Nstride0;
273                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
274                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride2+Nstride1           ;
275                  }
276                }
277              }
278              faceNECount+=local_NE1*local_NE0;
279         }         }
280           totalNECount+=NE1*NE0;
281       }       }
282       totalNECount+=NE1*NE2;       if (!periodic[0]  && (local_NE0>0)) {
283       faceNECount+=NE1*NE2;          /* **  elements on boundary 001 (x1=0): */
284    }      
285    if (!periodic[1]) {          if (e_offset0 == 0) {
286       /* **  elements on boundary 010 (x2=0): */             #pragma omp parallel for private(i1,i2,k,node0)
287                 for (i2=0;i2<local_NE2;i2++) {
288       #pragma omp parallel for private(i0,i2,k,node0)               for (i1=0;i1<local_NE1;i1++) {
289       for (i2=0;i2<NE2;i2++) {        
290         for (i0=0;i0<NE0;i0++) {                 k=i1+local_NE1*i2+faceNECount;
291           k=i0+NE0*i2+faceNECount;                 node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
292           node0=i0+N0*N1*i2;                 out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
293                   out->FaceElements->Tag[k]=1;
294                   out->FaceElements->Owner[k]=myRank;
295          
296                   if  (useElementsOnFace) {
297                      out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
298                      out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
299                      out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
300                      out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+Nstride1                      ;
301                      out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+Nstride0                      ;
302                      out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+Nstride2+Nstride0           ;
303                      out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+Nstride2+Nstride1+Nstride0;
304                      out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+Nstride1+Nstride0           ;
305                   } else {
306                      out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
307                      out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+Nstride2                      ;
308                      out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+Nstride2+Nstride1           ;
309                      out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           Nstride1           ;
310                   }
311                 }
312               }
313               faceNECount+=local_NE1*local_NE2;
314            }
315            totalNECount+=NE1*NE2;
316            /* **  elements on boundary 002 (x1=1): */
317            if (local_NE0+e_offset0 == NE0) {
318               #pragma omp parallel for private(i1,i2,k,node0)
319               for (i2=0;i2<local_NE2;i2++) {
320                 for (i1=0;i1<local_NE1;i1++) {
321                   k=i1+local_NE1*i2+faceNECount;
322                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
323                   out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
324                   out->FaceElements->Tag[k]=2;
325                   out->FaceElements->Owner[k]=myRank;
326        
327           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 if  (useElementsOnFace) {
328           out->FaceElements->Tag[k]=10;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      Nstride0;
329           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+16;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
330                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
331           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
332              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;        
333              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
334              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           Nstride1           ;
335              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1           ;
336              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2                      ;
337              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0+1;        
338              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N1*N0+N0+1;                 } else {
339              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0+N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +Nstride0;
340           } else {                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           Nstride1+Nstride0;
341              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
342              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2+           Nstride0;
343              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N1*N0+1;                 }
344              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N1*N0;          
345                 }
346               }
347               faceNECount+=local_NE1*local_NE2;
348           }           }
349         }           totalNECount+=NE1*NE2;
350       }       }
351       totalNECount+=NE0*NE2;       if (!periodic[1] && (local_NE1>0)) {
352       faceNECount+=NE0*NE2;          /* **  elements on boundary 010 (x2=0): */
353              if (e_offset1 == 0) {
354       /* **  elements on boundary 020 (x2=1): */             #pragma omp parallel for private(i0,i2,k,node0)
355                 for (i2=0;i2<local_NE2;i2++) {
356       #pragma omp parallel for private(i0,i2,k,node0)               for (i0=0;i0<local_NE0;i0++) {
357       for (i2=0;i2<NE2;i2++) {                 k=i0+local_NE0*i2+faceNECount;
358         for (i0=0;i0<NE0;i0++) {                 node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
359           k=i0+NE0*i2+faceNECount;          
360           node0=i0+(NE1-1)*N0+N0*N1*i2;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
361                   out->FaceElements->Tag[k]=10;
362                   out->FaceElements->Owner[k]=myRank;
363                   if  (useElementsOnFace) {
364                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
365                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
366                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2           +Nstride0;
367                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
368          
369                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           Nstride1           ;
370                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride1+           Nstride0;
371                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
372                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+Nstride2+Nstride1           ;
373                   } else {
374                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
375                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      Nstride0;
376                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+           Nstride0;
377                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride2                      ;
378                   }
379                 }
380               }
381               faceNECount+=local_NE0*local_NE2;
382            }
383            totalNECount+=NE0*NE2;
384            /* **  elements on boundary 020 (x2=1): */
385            if (local_NE1+e_offset1 == NE1) {
386               #pragma omp parallel for private(i0,i2,k,node0)
387               for (i2=0;i2<local_NE2;i2++) {
388                 for (i0=0;i0<local_NE0;i0++) {
389                   k=i0+local_NE0*i2+faceNECount;
390                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
391        
392           out->FaceElements->Tag[k]=20;                 out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
393           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 out->FaceElements->Tag[k]=20;
394           out->FaceElements->Color[k]=(i0%2)+2*(i2%2)+20;                 out->FaceElements->Owner[k]=myRank;
395          
396           if  (useElementsOnFace) {                 if  (useElementsOnFace) {
397              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
398              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
399              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
400              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+Nstride1+Nstride0           ;
401              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;        
402              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
403              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+1;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+Nstride2                      ;
404              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+Nstride2+           Nstride0;
405           } else {                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      Nstride0;
406              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+N0;                 } else {
407              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+N0*N1+N0;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           Nstride1           ;
408              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+N0*N1+N0+1;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+Nstride2+Nstride1           ;
409              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+N0+1;                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+Nstride2+Nstride1+Nstride0;
410           }                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           Nstride1+Nstride0;
411                   }
412         }               }
413               }
414               faceNECount+=local_NE0*local_NE2;
415            }
416            totalNECount+=NE0*NE2;
417       }       }
      totalNECount+=NE0*NE2;  
      faceNECount+=NE0*NE2;  
418    }    }
419    out->FaceElements->numColors=24;    if (Finley_noError()) {
420           /* add tag names */
421    /*  face elements done: */       Finley_Mesh_addTagMap(out,"top", 200);
422           Finley_Mesh_addTagMap(out,"bottom", 100);
423    /*   condense the nodes: */       Finley_Mesh_addTagMap(out,"left", 1);
424           Finley_Mesh_addTagMap(out,"right", 2);
425    Finley_Mesh_resolveNodeIds(out);       Finley_Mesh_addTagMap(out,"front", 10);
426         Finley_Mesh_addTagMap(out,"back", 20);
427      }
428    /* prepare mesh for further calculatuions:*/    /* prepare mesh for further calculatuions:*/
429    Finley_Mesh_prepare(out) ;    if (Finley_noError()) {
430             Finley_Mesh_resolveNodeIds(out);
431      }
432      if (Finley_noError()) {
433             Finley_Mesh_prepare(out, optimize);
434      }
435    
436    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    if (!Finley_noError()) {
437          Finley_Mesh_free(out);
438      }
439        /* free up memory */
440      Finley_ReferenceElementSet_dealloc(refPoints);
441      Finley_ReferenceElementSet_dealloc(refContactElements);
442      Finley_ReferenceElementSet_dealloc(refFaceElements);
443      Finley_ReferenceElementSet_dealloc(refElements);
444      Paso_MPIInfo_free( mpi_info );  
445    
   if (Finley_ErrorCode!=NO_ERROR) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
446    return out;    return out;
447  }  }
   
 /*  
 * $Log$  
 * Revision 1.1  2004/10/26 06:53:57  jgs  
 * Initial revision  
 *  
 * Revision 1.1.1.1  2004/06/24 04:00:40  johng  
 * Initial version of eys using boost-python.  
 *  
 *  
 */  
   

Legend:
Removed from v.82  
changed lines
  Added in v.2881

  ViewVC Help
Powered by ViewVC 1.1.26