/[escript]/branches/domexper/dudley/src/Mesh_hex20.c
ViewVC logotype

Diff of /branches/domexper/dudley/src/Mesh_hex20.c

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

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

Legend:
Removed from v.149  
changed lines
  Added in v.2748

  ViewVC Help
Powered by ViewVC 1.1.26