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

Legend:
Removed from v.155  
changed lines
  Added in v.1811

  ViewVC Help
Powered by ViewVC 1.1.26