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

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

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

revision 616 by elspeth, Wed Mar 22 02:46:56 2006 UTC revision 2722 by gross, Fri Oct 16 06:45:01 2009 UTC
# Line 1  Line 1 
1  /*  
2   ************************************************************  /*******************************************************
3   *          Copyright 2006 by ACcESS MNRF                   *  *
4   *                                                          *  * Copyright (c) 2003-2009 by University of Queensland
5   *              http://www.access.edu.au                    *  * Earth Systems Science Computational Center (ESSCC)
6   *       Primary Business: Queensland, Australia            *  * http://www.uq.edu.au/esscc
7   *  Licensed under the Open Software License version 3.0    *  *
8   *     http://www.opensource.org/licenses/osl-3.0.php       *  * 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    
# Line 21  Line 23 
23    
24  /**************************************************************/  /**************************************************************/
25    
 /*  Author: gross@access.edu.au */  
 /*  Version: $Id$  
   
 /**************************************************************/  
   
26  #include "RectangularMesh.h"  #include "RectangularMesh.h"
27    
28  /**************************************************************/  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,
29                                              double* Length,
30  Finley_Mesh* Finley_RectangularMesh_Hex20(dim_t* numElements,double* Length,bool_t* periodic,index_t order,bool_t useElementsOnFace) {                                            bool_t* periodic,
31    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,
32    index_t node0;                                            index_t reduced_order,
33                                              bool_t useElementsOnFace,
34                                              bool_t useFullElementOrder,
35                                              bool_t useMacroElements,
36                                              bool_t optimize)
37    {
38      #define N_PER_E 2
39      #define DIM 3
40      dim_t N0,N1,N2,NE0,NE1,NE2,i0,i1,i2,k,Nstride0,Nstride1,Nstride2, local_NE0, local_NE1, local_NE2;
41      dim_t totalNECount,faceNECount,NDOF0,NDOF1,NDOF2,NFaceElements, local_N0, local_N1, local_N2, NN;
42      index_t node0, myRank, e_offset0, e_offset1, e_offset2, offset0, offset1, offset2, global_i0, global_i1, global_i2;
43    Finley_Mesh* out;    Finley_Mesh* out;
44      Paso_MPIInfo *mpi_info = NULL;
45    char name[50];    char name[50];
46      bool_t generateAllNodes= useFullElementOrder || useMacroElements;
47      #ifdef Finley_TRACE
48    double time0=Finley_timer();    double time0=Finley_timer();
49      #endif
50    
51      /* get MPI information */
52      mpi_info = Paso_MPIInfo_alloc( MPI_COMM_WORLD );
53      if (! Finley_noError()) {
54            return NULL;
55      }
56      myRank=mpi_info->rank;
57    
58      /* set up the global dimensions of the mesh */
59    
60    NE0=MAX(1,numElements[0]);    NE0=MAX(1,numElements[0]);
61    NE1=MAX(1,numElements[1]);    NE1=MAX(1,numElements[1]);
62    NE2=MAX(1,numElements[2]);    NE2=MAX(1,numElements[2]);
63    N0=2*NE0+1;    N0=N_PER_E*NE0+1;
64    N1=2*NE1+1;    N1=N_PER_E*NE1+1;
65    N2=2*NE2+1;    N2=N_PER_E*NE2+1;
66    
67    if (N0<=MIN(N1,N2)) {    /*  allocate mesh: */  
68       if (N1 <= N2) {    sprintf(name,"Brick %d x %d x %d mesh",N0,N1,N2);
69          M0=1;    out=Finley_Mesh_alloc(name,DIM,order, reduced_order, mpi_info);
70          M1=N0;    if (! Finley_noError()) {
71          M2=N0*N1;        Paso_MPIInfo_free( mpi_info );
72          return NULL;
73      }
74    
75      if (generateAllNodes) {
76         /* Finley_setError(SYSTEM_ERROR,"full element order for Hex elements is not supported yet."); */
77         if (useMacroElements) {
78              Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex27Macro,
79                                                     out->order,
80                                                     out->reduced_order,
81                                                     mpi_info));
82       } else {       } else {
83          M0=1;            Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex27,
84          M2=N0;                                                   out->order,
85          M1=N0*N2;                                                   out->reduced_order,
86                                                     mpi_info));
87       }       }
88    } else if (N1<=MIN(N2,N0)) {       if (useElementsOnFace) {
89       if (N2 <= N0) {           Finley_setError(SYSTEM_ERROR,"rich elements for Hex27 elements is not supported yet.");
         M1=1;  
         M2=N1;  
         M0=N2*N1;  
90       } else {       } else {
91          M1=1;           if (useMacroElements) {
92          M0=N1;               Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec9Macro,
93          M2=N1*N0;                                                          out->order,
94                                                            out->reduced_order,
95                                                            mpi_info));
96             } else {
97                 Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec9,
98                                                            out->order,
99                                                            out->reduced_order,
100                                                            mpi_info));
101             }
102             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec9_Contact,
103                                                           out->order,
104                                                           out->reduced_order,
105                                                           mpi_info));
106       }       }
107    } else {  
108       if (N0 <= N1) {    } else  {
109          M2=1;       Finley_Mesh_setElements(out,Finley_ElementFile_alloc(Hex20,out->order,out->reduced_order,mpi_info));
110          M0=N2;       if (useElementsOnFace) {
111          M1=N2*N0;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Hex20Face,
112                                                                      out->order,
113                                                                      out->reduced_order,
114                                                                      mpi_info));
115             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Hex20Face_Contact,
116                                                                        out->order,
117                                                                        out->reduced_order,
118                                                                        mpi_info));
119       } else {       } else {
120          M2=1;           Finley_Mesh_setFaceElements(out,Finley_ElementFile_alloc(Rec8,
121          M1=N2;                                                                    out->order,
122          M0=N1*N2;                                                                    out->reduced_order,
123                                                                      mpi_info));
124             Finley_Mesh_setContactElements(out,Finley_ElementFile_alloc(Rec8_Contact,
125                                                                         out->order,
126                                                                         out->reduced_order,
127                                                                         mpi_info));
128       }       }
129    }    }
130      Finley_Mesh_setPoints(out,Finley_ElementFile_alloc(Point1,
131    NFaceElements=0;                                                   out->order,
132    if (!periodic[0]) {                                                   out->reduced_order,
133        NDOF0=N0;                                                   mpi_info));
134        NFaceElements+=2*NE1*NE2;    if (! Finley_noError()) {
135    } else {        Paso_MPIInfo_free( mpi_info );
136        NDOF0=N0-1;        Finley_Mesh_free(out);
137          return NULL;
138    }    }
139    if (!periodic[1]) {  
140        NDOF1=N1;    /* work out the largest dimension */
141        NFaceElements+=2*NE0*NE2;    if (N2==MAX3(N0,N1,N2)) {
142         Nstride0=1;
143         Nstride1=N0;
144         Nstride2=N0*N1;
145         local_NE0=NE0;
146         e_offset0=0;
147         local_NE1=NE1;
148         e_offset1=0;
149         Paso_MPIInfo_Split(mpi_info,NE2,&local_NE2,&e_offset2);
150      } else if (N1==MAX3(N0,N1,N2)) {
151         Nstride0=N2;
152         Nstride1=N0*N2;
153         Nstride2=1;
154         local_NE0=NE0;
155         e_offset0=0;
156         Paso_MPIInfo_Split(mpi_info,NE1,&local_NE1,&e_offset1);
157         local_NE2=NE2;
158         e_offset2=0;
159    } else {    } else {
160        NDOF1=N1-1;       Nstride0=N1*N2;
161    }       Nstride1=1;
162    if (!periodic[2]) {       Nstride2=N1;
163        NDOF2=N2;       Paso_MPIInfo_Split(mpi_info,NE0,&local_NE0,&e_offset0);
164        NFaceElements+=2*NE0*NE1;       local_NE1=NE1;
165         e_offset1=0;
166         local_NE2=NE2;
167         e_offset2=0;
168      }
169      offset0=e_offset0*N_PER_E;
170      offset1=e_offset1*N_PER_E;
171      offset2=e_offset2*N_PER_E;
172      local_N0=local_NE0>0 ? local_NE0*N_PER_E+1 : 0;
173      local_N1=local_NE1>0 ? local_NE1*N_PER_E+1 : 0;
174      local_N2=local_NE2>0 ? local_NE2*N_PER_E+1 : 0;
175    
176      /* get the number of surface elements */
177    
178      NFaceElements=0;
179      if (!periodic[2] && (local_NE2>0) ) {
180        NDOF2=N2;
181        if (offset2==0) NFaceElements+=local_NE1*local_NE0;
182        if (local_NE2+e_offset2 == NE2) NFaceElements+=local_NE1*local_NE0;
183    } else {    } else {
184        NDOF2=N2-1;        NDOF2=N2-1;
185    }    }
186      
187    /*  allocate mesh: */    if (!periodic[0] && (local_NE0>0) ) {
188           NDOF0=N0;
189    sprintf(name,"Rectangular %d x %d x %d mesh",N0,N1,N2);       if (e_offset0 == 0) NFaceElements+=local_NE1*local_NE2;
190    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(Hex20,out->order);  
   if (useElementsOnFace) {  
      out->FaceElements=Finley_ElementFile_alloc(Hex20Face,out->order);  
      out->ContactElements=Finley_ElementFile_alloc(Hex20Face_Contact,out->order);  
191    } else {    } else {
192       out->FaceElements=Finley_ElementFile_alloc(Rec8,out->order);        NDOF0=N0-1;
      out->ContactElements=Finley_ElementFile_alloc(Rec8_Contact,out->order);  
193    }    }
194    out->Points=Finley_ElementFile_alloc(Point1,out->order);    if (!periodic[1] && (local_NE1>0) ) {
195    if (! Finley_noError()) {       NDOF1=N1;
196        Finley_Mesh_dealloc(out);       if (e_offset1 == 0) NFaceElements+=local_NE0*local_NE2;
197        return NULL;       if (local_NE1+e_offset1 == NE1) NFaceElements+=local_NE0*local_NE2;
198      } else {
199          NDOF1=N1-1;
200    }    }
   
     
201    /*  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;  
   }  
202    
203    #pragma omp parallel for private(i0,i1,i2,k)    Finley_NodeFile_allocTable(out->Nodes,local_N0*local_N1*local_N2);
204    for (i2=0;i2<N2;i2++) {    Finley_ElementFile_allocTable(out->Elements,local_NE0*local_NE1*local_NE2);
205      for (i1=0;i1<N1;i1++) {    Finley_ElementFile_allocTable(out->FaceElements,NFaceElements);
       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);  
       }  
     }  
   }  
206    
207        if (Finley_noError()) {
208    /* tags for the faces: */       /* create nodes */
   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;  
       }  
     }  
   }  
     
   /*   set the elements: */  
     
   #pragma omp parallel for private(i0,i1,i2,k,node0)  
   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=2*i0+2*i1*N0+2*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,20)]=node0;  
         out->Elements->Nodes[INDEX2(1,k,20)]=node0+2;  
         out->Elements->Nodes[INDEX2(2,k,20)]=node0+2*N0+2;  
         out->Elements->Nodes[INDEX2(3,k,20)]=node0+2*N0;  
         out->Elements->Nodes[INDEX2(4,k,20)]=node0+2*N0*N1;  
         out->Elements->Nodes[INDEX2(5,k,20)]=node0+2*N0*N1+2;  
         out->Elements->Nodes[INDEX2(6,k,20)]=node0+2*N0*N1+2*N0+2;  
         out->Elements->Nodes[INDEX2(7,k,20)]=node0+2*N0*N1+2*N0;  
         out->Elements->Nodes[INDEX2(8,k,20)]=node0+1;  
         out->Elements->Nodes[INDEX2(9,k,20)]=node0+N0+2;  
         out->Elements->Nodes[INDEX2(10,k,20)]=node0+2*N0+1;  
         out->Elements->Nodes[INDEX2(11,k,20)]=node0+N0;  
         out->Elements->Nodes[INDEX2(12,k,20)]=node0+N0*N1;  
         out->Elements->Nodes[INDEX2(13,k,20)]=node0+N0*N1+2;  
         out->Elements->Nodes[INDEX2(14,k,20)]=node0+N0*N1+2*N0+2;  
         out->Elements->Nodes[INDEX2(15,k,20)]=node0+N0*N1+2*N0;  
         out->Elements->Nodes[INDEX2(16,k,20)]=node0+2*N0*N1+1;  
         out->Elements->Nodes[INDEX2(17,k,20)]=node0+2*N0*N1+N0+2;  
         out->Elements->Nodes[INDEX2(18,k,20)]=node0+2*N0*N1+2*N0+1;  
         out->Elements->Nodes[INDEX2(19,k,20)]=node0+2*N0*N1+N0;  
       }  
     }  
   }  
   out->Elements->minColor=0;  
   out->Elements->maxColor=COLOR_MOD(0)+3*COLOR_MOD(0)+9*COLOR_MOD(0);  
     
   /*   face elements: */  
     
   if  (useElementsOnFace) {  
      NUMNODES=20;  
   } else {  
      NUMNODES=8;  
   }  
   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=2*i0+2*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+2*N0;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0;  
            out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+1;  
            out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;  
            out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0*N1+N0;  
            out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+2*N0*N1+1;  
         } else {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(7,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=2*i0+2*i1*N0+2*N0*N1*(NE2-1);  
     
         out->FaceElements->Id[k]=i0+NE0*i1+totalNECount;  
         out->FaceElements->Tag[k]=200;  
         out->FaceElements->Color[k]=(i0%2)+2*(i1%2)+4;  
   
         if  (useElementsOnFace) {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+2*N0*N1+1;  
            out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0*N1+N0;  
   
            out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0*N1;  
            out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0*N1+2*N0;  
   
            out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+1;  
            out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0+2;  
            out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0;  
   
         } else {  
            out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0*N1;  
            out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2;  
            out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
            out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
            out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0*N1+1;  
            out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
            out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
            out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1+N0;  
         }  
       }  
     }  
     totalNECount+=NE1*NE0;  
     faceNECount+=NE1*NE0;  
   }  
   if (!periodic[0]) {  
      /* **  elements on boundary 001 (x1=0): */  
     
      #pragma omp parallel for private(i1,i2,k,node0)  
      for (i2=0;i2<NE2;i2++) {  
        for (i1=0;i1<NE1;i1++) {  
          k=i1+NE1*i2+faceNECount;  
          node0=2*i1*N0+2*N0*N1*i2;  
209        
210           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;       #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
211           out->FaceElements->Tag[k]=1;       for (i2=0;i2<local_N2;i2++) {
212           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+8;         for (i1=0;i1<local_N1;i1++) {
213             for (i0=0;i0<local_N0;i0++) {
214           if  (useElementsOnFace) {             k=i0+local_N0*i1+local_N0*local_N1*i2;
215              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;             global_i0=i0+offset0;
216              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;             global_i1=i1+offset1;
217              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;             global_i2=i2+offset2;
218              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;             out->Nodes->Coordinates[INDEX2(0,k,DIM)]=DBLE(global_i0)/DBLE(N0-1)*Length[0];
219               out->Nodes->Coordinates[INDEX2(1,k,DIM)]=DBLE(global_i1)/DBLE(N1-1)*Length[1];
220              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2;             out->Nodes->Coordinates[INDEX2(2,k,DIM)]=DBLE(global_i2)/DBLE(N2-1)*Length[2];
221              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+2;             out->Nodes->Id[k]=Nstride0*global_i0+Nstride1*global_i1+Nstride2*global_i2;
222              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;             out->Nodes->Tag[k]=0;
223              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+2;             out->Nodes->globalDegreesOfFreedom[k]=Nstride0*(global_i0%NDOF0)
224                                                   +Nstride1*(global_i1%NDOF1)
225              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0*N1;                                                 +Nstride2*(global_i2%NDOF2);
             out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0;  
   
             out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;  
             out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+1;  
             out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;  
             out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0+1;  
   
             out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
             out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0+2;  
   
          } else {  
             out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;  
             out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1;  
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0*N1;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1+N0;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0;  
226           }           }
227         }         }
228       }       }
229       totalNECount+=NE1*NE2;       /*   set the elements: */
230       faceNECount+=NE1*NE2;       NN=out->Elements->numNodes;
231           #pragma omp parallel for private(i0,i1,i2,k,node0)
232       /* **  elements on boundary 002 (x1=1): */       for (i2=0;i2<local_NE2;i2++) {
233             for (i1=0;i1<local_NE1;i1++) {
234       #pragma omp parallel for private(i1,i2,k,node0)           for (i0=0;i0<local_NE0;i0++) {
235       for (i2=0;i2<NE2;i2++) {            
236         for (i1=0;i1<NE1;i1++) {             k=i0+local_NE0*i1+local_NE0*local_NE1*i2;        
237           k=i1+NE1*i2+faceNECount;             node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
          node0=2*(NE0-1)+2*i1*N0+2*N0*N1*i2 ;  
238        
239           out->FaceElements->Id[k]=i1+NE1*i2+totalNECount;             out->Elements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+NE0*NE1*(i2+e_offset2);
240           out->FaceElements->Tag[k]=2;             out->Elements->Tag[k]=0;
241           out->FaceElements->Color[k]=(i2%2)+2*(i1%2)+12;             out->Elements->Owner[k]=myRank;
242    
243           if  (useElementsOnFace) {             out->Elements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
244              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;             out->Elements->Nodes[INDEX2(1,k,NN)] =node0+                      2*Nstride0;
245              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;             out->Elements->Nodes[INDEX2(2,k,NN)] =node0+           2*Nstride1+2*Nstride0;
246              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;             out->Elements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride1;
247              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;             out->Elements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2                      ;
248               out->Elements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2           +2*Nstride0;
249              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;             out->Elements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
250              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0;             out->Elements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
251              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2*N0;             out->Elements->Nodes[INDEX2(8,k,NN)] =node0+                      1*Nstride0;
252              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0*N1;             out->Elements->Nodes[INDEX2(9,k,NN)] =node0+           1*Nstride1+2*Nstride0;
253               out->Elements->Nodes[INDEX2(10,k,NN)]=node0+           2*Nstride1+1*Nstride0;
254              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N0+2;             out->Elements->Nodes[INDEX2(11,k,NN)]=node0+           1*Nstride1           ;
255              out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2*N0+2;             out->Elements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2                      ;
256              out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N0*N1+N0+2;             out->Elements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2           +2*Nstride0;
257              out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N0*N1+2;             out->Elements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
258               out->Elements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
259              out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+1;             out->Elements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2           +1*Nstride0;
260              out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0+1;             out->Elements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
261              out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+2*N0+1;             out->Elements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
262              out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N0*N1+1;             out->Elements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
263               if (generateAllNodes) {
264              out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N0;                out->Elements->Nodes[INDEX2(20,k,NN)]=node0+           1*Nstride1+1*Nstride0;
265              out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0;                out->Elements->Nodes[INDEX2(21,k,NN)]=node0+1*Nstride2           +1*Nstride0;
266              out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N0*N1+N0;                out->Elements->Nodes[INDEX2(22,k,NN)]=node0+1*Nstride2+1*Nstride1+2*Nstride0;
267              out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N0*N1;                out->Elements->Nodes[INDEX2(23,k,NN)]=node0+1*Nstride2+2*Nstride1+1*Nstride0;
268                  out->Elements->Nodes[INDEX2(24,k,NN)]=node0+1*Nstride2+1*Nstride1           ;
269           } else {                out->Elements->Nodes[INDEX2(25,k,NN)]=node0+2*Nstride2+1*Nstride1+1*Nstride0;
270              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2;                out->Elements->Nodes[INDEX2(26,k,NN)]=node0+1*Nstride2+1*Nstride1+1*Nstride0;        
271              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0+2;             }
             out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0*N1+2;  
             out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N0+2;  
             out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2*N0+2;  
             out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+N0+2;  
             out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N0*N1+2;  
272           }           }
     
273         }         }
274       }       }
275       totalNECount+=NE1*NE2;       /* face elements */
276       faceNECount+=NE1*NE2;       NN=out->FaceElements->numNodes;
277    }       totalNECount=NE0*NE1*NE2;
278    if (!periodic[1]) {       faceNECount=0;
279       /* **  elements on boundary 010 (x2=0): */       /*   these are the quadrilateral elements on boundary 1 (x3=0): */
280           if (!periodic[2] && (local_NE2>0)) {
281       #pragma omp parallel for private(i0,i2,k,node0)         /* **  elements on boundary 100 (x3=0): */
282       for (i2=0;i2<NE2;i2++) {         if (offset2==0) {
283         for (i0=0;i0<NE0;i0++) {            #pragma omp parallel for private(i0,i1,k,node0)
284           k=i0+NE0*i2+faceNECount;            for (i1=0;i1<local_NE1;i1++) {
285           node0=2*i0+2*N0*N1*i2;              for (i0=0;i0<local_NE0;i0++) {
286                
287           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                k=i0+local_NE0*i1+faceNECount;
288           out->FaceElements->Tag[k]=10;                node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1);
289           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+16;      
290                  out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
291           if  (useElementsOnFace) {                out->FaceElements->Tag[k]=100;
292              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                out->FaceElements->Owner[k]=myRank;
293              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;          
294              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;                if  (useElementsOnFace) {
295              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;                   out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
296                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0           +2*Nstride1           ;
297              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+2*N0;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0           +2*Nstride1+2*Nstride0;
298              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0+2;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride0           ;
299              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+2*N0+2;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2                      ;
300              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N1*N0+2*N0;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
301                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
302              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+1;                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2           +2*Nstride0;
303              out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+N0*N1+2;                   out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+           1*Nstride1           ;
304              out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+2*N1*N0+1;                   out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+           2*Nstride1+1*Nstride0;
305              out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+N1*N0;                   out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+           1*Nstride1+2*Nstride0;
306                     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+                      1*Nstride0;
307              out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;                   out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2                      ;
308              out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+N0+2;                   out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
309              out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N1*N0+N0+2;                   out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
310              out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+2*N1*N0+N0;                   out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2           +2*Nstride0;
311                     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+2*Nstride2+1*Nstride1;
312              out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+2*N0+1;                   out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
313              out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+N0*N1+2*N0+2;                   out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
314              out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;                   out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+2*Nstride2           +1*Nstride0;
315              out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+N1*N0+2*N0;                } else {
316                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
317           } else {                   out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+           2*Nstride1           ;
318              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0;                   out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+           2*Nstride1+2*Nstride0;
319              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2;                   out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+                      2*Nstride0;
320              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N1*N0+2;                   out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+           1*Nstride1           ;
321              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N1*N0;                   out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+           2*Nstride1+1*Nstride0;
322              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+1;                   out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           1*Nstride1+2*Nstride0;
323              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+N0*N1+2;                   out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+                      1*Nstride0;
324              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N1*N0+1;                   if (generateAllNodes){
325              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+N1*N0;                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+           1*Nstride1+1*Nstride0;
326           }                   }
327                  }
328                }
329              }
330              faceNECount+=local_NE1*local_NE0;
331           }
332           totalNECount+=NE1*NE0;
333           /* **  elements on boundary 200 (x3=1) */
334           if (local_NE2+e_offset2 == NE2) {
335              #pragma omp parallel for private(i0,i1,k,node0)
336              for (i1=0;i1<local_NE1;i1++) {
337                for (i0=0;i0<local_NE0;i0++) {
338          
339                  k=i0+local_NE0*i1+faceNECount;
340                  node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(NE2-1);
341            
342                  out->FaceElements->Id[k]=(i0+e_offset0)+NE0*(i1+e_offset1)+totalNECount;
343                  out->FaceElements->Tag[k]=200;
344                  out->FaceElements->Owner[k]=myRank;
345                  if  (useElementsOnFace) {
346                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2                      ;
347                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2+           2*Nstride0;
348                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
349                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
350          
351                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0                                 ;
352                     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride0                      ;
353                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+           2*Nstride1+2*Nstride0;
354                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           2*Nstride1;
355          
356                     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+           1*Nstride0;
357                     out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
358                     out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
359                     out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
360          
361                     out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+1*Nstride2;
362                     out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+1*Nstride2           +2*Nstride0;
363                     out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
364                     out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
365          
366                     out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+                      1*Nstride0;
367                     out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+           1*Nstride1+2*Nstride0;
368                     out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+           2*Nstride1+1*Nstride0;
369                     out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+           1*Nstride1           ;
370          
371                  } else {
372                     out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0+2*Nstride2                      ;
373                     out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2           +2*Nstride0;
374                     out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
375                     out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
376                     out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride2           +1*Nstride0;
377                     out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1+2*Nstride0;
378                     out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+1*Nstride0;
379                     out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
380                     if (generateAllNodes){
381                     out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+2*Nstride2+1*Nstride1+1*Nstride0;
382                     }
383                  }
384                }
385              }
386              faceNECount+=local_NE1*local_NE0;
387         }         }
388           totalNECount+=NE1*NE0;
389       }       }
390       totalNECount+=NE0*NE2;       if (!periodic[0] && (local_NE0>0)) {
391       faceNECount+=NE0*NE2;          /* **  elements on boundary 001 (x1=0): */
392          
393       /* **  elements on boundary 020 (x2=1): */          if (e_offset0 == 0) {
394               #pragma omp parallel for private(i1,i2,k,node0)
395               for (i2=0;i2<local_NE2;i2++) {
396                 for (i1=0;i1<local_NE1;i1++) {
397          
398                   k=i1+local_NE1*i2+faceNECount;
399                   node0=Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
400                   out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
401                   out->FaceElements->Tag[k]=1;
402                   out->FaceElements->Owner[k]=myRank;
403          
404                   if  (useElementsOnFace) {
405                      out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
406                      out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2                      ;
407                      out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
408                      out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+2*Nstride1                      ;
409          
410                      out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+2*Nstride0                      ;
411                      out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+2*Nstride0           ;
412                      out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+2*Nstride2+2*Nstride1+2*Nstride0;
413                      out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+2*Nstride1+2*Nstride0           ;
414          
415                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2                      ;
416                      out->FaceElements->Nodes[INDEX2(9,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
417                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
418                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+           1*Nstride1           ;
419          
420                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+                      1*Nstride0;
421                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2           +1*Nstride0;
422                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
423                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride1+           1*Nstride0;
424          
425                      out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2+           2*Nstride0;
426                      out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
427                      out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
428                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride1+           2*Nstride0;
429          
430                   } else {
431                      out->FaceElements->Nodes[INDEX2(0,k,NN)] =node0                                 ;
432                      out->FaceElements->Nodes[INDEX2(1,k,NN)] =node0+2*Nstride2                      ;
433                      out->FaceElements->Nodes[INDEX2(2,k,NN)] =node0+2*Nstride2+2*Nstride1           ;
434                      out->FaceElements->Nodes[INDEX2(3,k,NN)] =node0+           2*Nstride1           ;
435                      out->FaceElements->Nodes[INDEX2(4,k,NN)] =node0+1*Nstride2                      ;
436                      out->FaceElements->Nodes[INDEX2(5,k,NN)] =node0+2*Nstride2+1*Nstride1           ;
437                      out->FaceElements->Nodes[INDEX2(6,k,NN)] =node0+1*Nstride2+2*Nstride1           ;
438                      out->FaceElements->Nodes[INDEX2(7,k,NN)] =node0+           1*Nstride1           ;
439                     if (generateAllNodes){
440                        out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1           ;
441                     }
442                   }
443                 }
444               }
445               faceNECount+=local_NE1*local_NE2;
446            }
447            totalNECount+=NE1*NE2;
448            
449       #pragma omp parallel for private(i0,i2,k,node0)          /* **  elements on boundary 002 (x1=1): */
450       for (i2=0;i2<NE2;i2++) {          if (local_NE0+e_offset0 == NE0) {
451         for (i0=0;i0<NE0;i0++) {             #pragma omp parallel for private(i1,i2,k,node0)
452           k=i0+NE0*i2+faceNECount;             for (i2=0;i2<local_NE2;i2++) {
453           node0=2*i0+2*(NE1-1)*N0+2*N0*N1*i2;               for (i1=0;i1<local_NE1;i1++) {
454                   k=i1+local_NE1*i2+faceNECount;
455                   node0=Nstride0*N_PER_E*(NE0-1)+Nstride1*N_PER_E*(i1+e_offset1)+Nstride2*N_PER_E*(i2+e_offset2);
456                   out->FaceElements->Id[k]=(i1+e_offset1)+NE1*(i2+e_offset2)+totalNECount;
457                   out->FaceElements->Tag[k]=2;
458                   out->FaceElements->Owner[k]=myRank;
459        
460           out->FaceElements->Id[k]=i2+NE2*i0+totalNECount;                 if  (useElementsOnFace) {
461           out->FaceElements->Tag[k]=20;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+                      2*Nstride0;
462           out->FaceElements->Color[k]=(i2%2)+2*(i0%2)+20;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           2*Nstride1+2*Nstride0;
463                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
464           if  (useElementsOnFace) {                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+           2*Nstride0;
465              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;        
466              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
467              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+           2*Nstride1           ;
468              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
469                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2                      ;
470              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0;        
471              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N0*N1;                    out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+           1*Nstride1+2*Nstride0;
472              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+2*N0*N1+2;                    out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
473              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2;                    out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
474                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2+           2*Nstride0;
475              out->FaceElements->Nodes[INDEX2(8,k,NUMNODES)]=node0+N1*N0+2*N0;        
476              out->FaceElements->Nodes[INDEX2(9,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;                    out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+                      1*Nstride0;
477              out->FaceElements->Nodes[INDEX2(10,k,NUMNODES)]=node0+N0*N1+2*N0+2;                    out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+           2*Nstride1+1*Nstride0;
478              out->FaceElements->Nodes[INDEX2(11,k,NUMNODES)]=node0+2*N0+1;                    out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
479                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+           1*Nstride0;
480              out->FaceElements->Nodes[INDEX2(12,k,NUMNODES)]=node0+N0;        
481              out->FaceElements->Nodes[INDEX2(13,k,NUMNODES)]=node0+2*N0*N1+N0;                    out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+           1*Nstride1           ;
482              out->FaceElements->Nodes[INDEX2(14,k,NUMNODES)]=node0+2*N0*N1+N0+2;                    out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
483              out->FaceElements->Nodes[INDEX2(15,k,NUMNODES)]=node0+N0+2;                    out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
484                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2                      ;
485              out->FaceElements->Nodes[INDEX2(16,k,NUMNODES)]=node0+N1*N0;        
486              out->FaceElements->Nodes[INDEX2(17,k,NUMNODES)]=node0+2*N1*N0+1;                 } else {
487              out->FaceElements->Nodes[INDEX2(18,k,NUMNODES)]=node0+N0*N1+2;                    out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                      +2*Nstride0;
488              out->FaceElements->Nodes[INDEX2(19,k,NUMNODES)]=node0+1;                    out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+           2*Nstride1+2*Nstride0;
489           } else {                    out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
490              out->FaceElements->Nodes[INDEX2(0,k,NUMNODES)]=node0+2*N0;                    out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2+           2*Nstride0;
491              out->FaceElements->Nodes[INDEX2(1,k,NUMNODES)]=node0+2*N0*N1+2*N0;                    out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           1*Nstride1+2*Nstride0;
492              out->FaceElements->Nodes[INDEX2(2,k,NUMNODES)]=node0+2*N0*N1+2*N0+2;                    out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
493              out->FaceElements->Nodes[INDEX2(3,k,NUMNODES)]=node0+2*N0+2;                    out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
494              out->FaceElements->Nodes[INDEX2(4,k,NUMNODES)]=node0+N1*N0+2*N0;                    out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2           +2*Nstride0;
495              out->FaceElements->Nodes[INDEX2(5,k,NUMNODES)]=node0+2*N1*N0+2*N0+1;                   if (generateAllNodes){
496              out->FaceElements->Nodes[INDEX2(6,k,NUMNODES)]=node0+N0*N1+2*N0+2;                      out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+1*Nstride1+2*Nstride0;
497              out->FaceElements->Nodes[INDEX2(7,k,NUMNODES)]=node0+2*N0+1;                   }
498                   }
499            
500                 }
501               }
502               faceNECount+=local_NE1*local_NE2;
503           }           }
504         }           totalNECount+=NE1*NE2;
505         }
506         if (!periodic[1] && (local_NE1>0)) {
507            /* **  elements on boundary 010 (x2=0): */
508            if (e_offset1 == 0) {
509               #pragma omp parallel for private(i0,i2,k,node0)
510               for (i2=0;i2<local_NE2;i2++) {
511                 for (i0=0;i0<local_NE0;i0++) {
512                   k=i0+local_NE0*i2+faceNECount;
513                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride2*N_PER_E*(i2+e_offset2);
514            
515                   out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(e_offset0+i0)+totalNECount;
516                   out->FaceElements->Tag[k]=10;
517                   out->FaceElements->Owner[k]=myRank;
518                   if  (useElementsOnFace) {
519                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
520                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      2*Nstride0;
521                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2           +2*Nstride0;
522                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2                      ;
523          
524                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+           2*Nstride1           ;
525                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride1+           2*Nstride0;
526                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
527                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
528          
529                      out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+                      1*Nstride0;
530                      out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+1*Nstride2+           2*Nstride0;
531                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+2*Nstride2+           1*Nstride0;
532                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+1*Nstride2                      ;
533          
534                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+           1*Nstride1           ;
535                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+           1*Nstride1+2*Nstride0;
536                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
537                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
538      
539                      out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+           2*Nstride1+1*Nstride0;
540                      out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
541                      out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
542                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
543          
544                   } else {
545                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0                                 ;
546                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+                      2*Nstride0;
547                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+           2*Nstride0;
548                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride2                      ;
549                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+                      1*Nstride0;
550                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+1*Nstride2+           2*Nstride0;
551                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+           1*Nstride0;
552                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+1*Nstride2                      ;
553                     if (generateAllNodes){
554                        out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+         1*Nstride0;
555                     }
556                   }
557                 }
558               }
559               faceNECount+=local_NE0*local_NE2;
560            }
561            totalNECount+=NE0*NE2;
562            /* **  elements on boundary 020 (x2=1): */
563            if (local_NE1+e_offset1 == NE1) {
564               #pragma omp parallel for private(i0,i2,k,node0)
565               for (i2=0;i2<local_NE2;i2++) {
566                 for (i0=0;i0<local_NE0;i0++) {
567                   k=i0+local_NE0*i2+faceNECount;
568                   node0=Nstride0*N_PER_E*(i0+e_offset0)+Nstride1*N_PER_E*(NE1-1)+Nstride2*N_PER_E*(i2+e_offset2);
569      
570                   out->FaceElements->Id[k]=(i2+e_offset2)+NE2*(i0+e_offset0)+totalNECount;
571                   out->FaceElements->Tag[k]=20;
572                   out->FaceElements->Owner[k]=myRank;
573          
574                   if  (useElementsOnFace) {
575                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           2*Nstride1           ;
576                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
577                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
578                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+2*Nstride1+2*Nstride0           ;
579          
580                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0                                 ;
581                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2                      ;
582                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+2*Nstride2+           2*Nstride0;
583                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+                      2*Nstride0;
584          
585                      out->FaceElements->Nodes[INDEX2(8,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
586                      out->FaceElements->Nodes[INDEX2(9,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
587                      out->FaceElements->Nodes[INDEX2(10,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
588                      out->FaceElements->Nodes[INDEX2(11,k,NN)]=node0+           2*Nstride1+1*Nstride0;
589          
590                      out->FaceElements->Nodes[INDEX2(12,k,NN)]=node0+           1*Nstride1           ;
591                      out->FaceElements->Nodes[INDEX2(13,k,NN)]=node0+2*Nstride2+1*Nstride1           ;
592                      out->FaceElements->Nodes[INDEX2(14,k,NN)]=node0+2*Nstride2+1*Nstride1+2*Nstride0;
593                      out->FaceElements->Nodes[INDEX2(15,k,NN)]=node0+           1*Nstride1+2*Nstride0;
594          
595                      out->FaceElements->Nodes[INDEX2(16,k,NN)]=node0+1*Nstride2                      ;
596                      out->FaceElements->Nodes[INDEX2(17,k,NN)]=node0+2*Nstride2           +1*Nstride0;
597                      out->FaceElements->Nodes[INDEX2(18,k,NN)]=node0+1*Nstride2           +2*Nstride0;
598                      out->FaceElements->Nodes[INDEX2(19,k,NN)]=node0+                      1*Nstride0;
599                   } else {
600                      out->FaceElements->Nodes[INDEX2(0,k,NN)]=node0+           2*Nstride1           ;
601                      out->FaceElements->Nodes[INDEX2(1,k,NN)]=node0+2*Nstride2+2*Nstride1           ;
602                      out->FaceElements->Nodes[INDEX2(2,k,NN)]=node0+2*Nstride2+2*Nstride1+2*Nstride0;
603                      out->FaceElements->Nodes[INDEX2(3,k,NN)]=node0+           2*Nstride1+2*Nstride0;
604                      out->FaceElements->Nodes[INDEX2(4,k,NN)]=node0+1*Nstride2+2*Nstride1           ;
605                      out->FaceElements->Nodes[INDEX2(5,k,NN)]=node0+2*Nstride2+2*Nstride1+1*Nstride0;
606                      out->FaceElements->Nodes[INDEX2(6,k,NN)]=node0+1*Nstride2+2*Nstride1+2*Nstride0;
607                      out->FaceElements->Nodes[INDEX2(7,k,NN)]=node0+           2*Nstride1+1*Nstride0;
608                     if (generateAllNodes){
609                        out->FaceElements->Nodes[INDEX2(8,k,NN)] =node0+1*Nstride2+2*Nstride1+1*Nstride0;
610                     }
611                   }
612                 }
613               }
614               faceNECount+=local_NE0*local_NE2;
615            }
616            totalNECount+=NE0*NE2;
617         }
618         /* add tag names */
619         Finley_Mesh_addTagMap(out,"top", 200);
620         Finley_Mesh_addTagMap(out,"bottom", 100);
621         Finley_Mesh_addTagMap(out,"left", 1);
622         Finley_Mesh_addTagMap(out,"right", 2);
623         Finley_Mesh_addTagMap(out,"front", 10);
624         Finley_Mesh_addTagMap(out,"back", 20);
625      
626         /* prepare mesh for further calculatuions:*/
627         if (Finley_noError()) {
628             Finley_Mesh_resolveNodeIds(out);
629         }
630         if (Finley_noError()) {
631             Finley_Mesh_prepare(out, optimize);
632       }       }
      totalNECount+=NE0*NE2;  
      faceNECount+=NE0*NE2;  
633    }    }
   out->FaceElements->minColor=0;  
   out->FaceElements->maxColor=24;  
   
   /*  face elements done: */  
     
   /*   condense the nodes: */  
     
   Finley_Mesh_resolveNodeIds(out);  
   
   /* prepare mesh for further calculatuions:*/  
   Finley_Mesh_prepare(out) ;  
634    
635      if (!Finley_noError()) {
636          Finley_Mesh_free(out);
637      }
638      /* free up memory */
639      Paso_MPIInfo_free( mpi_info );  
640    #ifdef Finley_TRACE    #ifdef Finley_TRACE
641    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);    printf("timing: mesh generation: %.4e sec\n",Finley_timer()-time0);
642    #endif    #endif
643    
   if (! Finley_noError()) {  
       Finley_Mesh_dealloc(out);  
       return NULL;  
   }  
644    return out;    return out;
645  }  }
   

Legend:
Removed from v.616  
changed lines
  Added in v.2722

  ViewVC Help
Powered by ViewVC 1.1.26