/[escript]/branches/trilinos_from_5897/dudley/src/Mesh_tet4.cpp
ViewVC logotype

Diff of /branches/trilinos_from_5897/dudley/src/Mesh_tet4.cpp

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

revision 6008 by caltinay, Mon Feb 22 06:59:27 2016 UTC revision 6009 by caltinay, Wed Mar 2 04:13:26 2016 UTC
# Line 14  Line 14 
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16    
17  /************************************************************************************/  /****************************************************************************/
18    
19  /*   Dudley: generates rectangular meshes */  /*   Dudley: generates rectangular meshes */
20    
# Line 22  Line 22 
22  /*   [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */  /*   [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
23  /*   integration scheme. */  /*   integration scheme. */
24    
25  /************************************************************************************/  /****************************************************************************/
   
 #define ESNEEDPYTHON  
 #include "esysUtils/first.h"  
26    
27  #include "TriangularMesh.h"  #include "TriangularMesh.h"
28    
29  /* Be careful reading this function. The X? and NStride? are 1,2,3 but the loop vars are 0,1,2 */  #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_))
30    
31    namespace dudley {
32    
33    // Be careful reading this function. The X? and NStride? are 1,2,3
34    // but the loop vars are 0,1,2
35  Dudley_Mesh *Dudley_TriangularMesh_Tet4(dim_t * numElements,  Dudley_Mesh *Dudley_TriangularMesh_Tet4(dim_t * numElements,
36                      double *Length, index_t order, index_t reduced_order, bool optimize, esysUtils::JMPI& mpi_info)                                          double *Length, index_t order, index_t reduced_order, bool optimize, escript::JMPI& mpi_info)
37  {  {
38  #define N_PER_E 1  #define N_PER_E 1
39  #define DIM 3  #define DIM 3
40      dim_t N0, N1, N2, NE0, NE1, NE2, i0, i1, i2, k, Nstride0 = 0, Nstride1 = 0, Nstride2 =      dim_t N0, N1, N2, NE0, NE1, NE2, i0, i1, i2, k, Nstride0 = 0, Nstride1 = 0, Nstride2 =
41      0, local_NE0, local_NE1, local_NE2, local_N0 = 0, local_N1 = 0, local_N2 = 0;          0, local_NE0, local_NE1, local_NE2, local_N0 = 0, local_N1 = 0, local_N2 = 0;
42      dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NDOF2 = 0, NFaceElements = 0, NN;      dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NDOF2 = 0, NFaceElements = 0, NN;
43      index_t node0, myRank, e_offset2, e_offset1, e_offset0 = 0, offset1 = 0, offset2 = 0, offset0 =      index_t node0, myRank, e_offset2, e_offset1, e_offset0 = 0, offset1 = 0, offset2 = 0, offset0 =
44      0, global_i0, global_i1, global_i2;          0, global_i0, global_i1, global_i2;
45      Dudley_Mesh *out;      Dudley_Mesh *out;
46      char name[50];      char name[50];
47  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
48      double time0 = Dudley_timer();      double time0 = Dudley_timer();
49  #endif  #endif
50    
51      const int LEFTTAG = 1;  /* boundary x1=0 */      const int LEFTTAG = 1;      /* boundary x1=0 */
52      const int RIGHTTAG = 2; /* boundary x1=1 */      const int RIGHTTAG = 2;     /* boundary x1=1 */
53      const int BOTTOMTAG = 100;  /* boundary x3=1 */      const int BOTTOMTAG = 100;  /* boundary x3=1 */
54      const int TOPTAG = 200; /* boundary x3=0 */      const int TOPTAG = 200;     /* boundary x3=0 */
55      const int FRONTTAG = 10;    /* boundary x2=0 */      const int FRONTTAG = 10;    /* boundary x2=0 */
56      const int BACKTAG = 20; /* boundary x2=1 */      const int BACKTAG = 20;     /* boundary x2=1 */
57    
58      /* get MPI information */      /* get MPI information */
59      myRank = mpi_info->rank;      myRank = mpi_info->rank;
60    
61      /* set up the global dimensions of the mesh */      /* set up the global dimensions of the mesh */
62        NE0 = std::max(1, numElements[0]);
63      NE0 = MAX(1, numElements[0]);      NE1 = std::max(1, numElements[1]);
64      NE1 = MAX(1, numElements[1]);      NE2 = std::max(1, numElements[2]);
     NE2 = MAX(1, numElements[2]);  
65      N0 = N_PER_E * NE0 + 1;      N0 = N_PER_E * NE0 + 1;
66      N1 = N_PER_E * NE1 + 1;      N1 = N_PER_E * NE1 + 1;
67      N2 = N_PER_E * NE2 + 1;      N2 = N_PER_E * NE2 + 1;
# Line 68  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 69  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
69      /*  allocate mesh: */      /*  allocate mesh: */
70      sprintf(name, "Triangular %d x %d x %d (x 5) mesh", N0, N1, N2);      sprintf(name, "Triangular %d x %d x %d (x 5) mesh", N0, N1, N2);
71      out = Dudley_Mesh_alloc(name, DIM, mpi_info);      out = Dudley_Mesh_alloc(name, DIM, mpi_info);
     if (!Dudley_noError())  
     {  
     return NULL;  
     }  
     if (Dudley_noError())  
     {  
72    
73      Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));      Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));
74      Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));      Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));
75      Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tet4, mpi_info));      Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tet4, mpi_info));
76    
77      /* work out the largest dimension */      /* work out the largest dimension */
78      if (N2 == MAX3(N0, N1, N2))      if (N2 == MAX3(N0, N1, N2)) {
79      {          Nstride0 = 1;
80          Nstride0 = 1;          Nstride1 = N0;
81          Nstride1 = N0;          Nstride2 = N0 * N1;
82          Nstride2 = N0 * N1;          local_NE0 = NE0;
83          local_NE0 = NE0;          e_offset0 = 0;
84          e_offset0 = 0;          local_NE1 = NE1;
85          local_NE1 = NE1;          e_offset1 = 0;
86          e_offset1 = 0;          mpi_info->split(NE2, &local_NE2, &e_offset2);
87          mpi_info->split(NE2, &local_NE2, &e_offset2);      } else if (N1 == MAX3(N0, N1, N2)) {
88      }          Nstride0 = N2;
89      else if (N1 == MAX3(N0, N1, N2))          Nstride1 = N0 * N2;
90      {          Nstride2 = 1;
91          Nstride0 = N2;          local_NE0 = NE0;
92          Nstride1 = N0 * N2;          e_offset0 = 0;
93          Nstride2 = 1;          mpi_info->split(NE1, &local_NE1, &e_offset1);
94          local_NE0 = NE0;          local_NE2 = NE2;
95          e_offset0 = 0;          e_offset2 = 0;
96          mpi_info->split(NE1, &local_NE1, &e_offset1);      } else {
97          local_NE2 = NE2;          Nstride0 = N1 * N2;
98          e_offset2 = 0;          Nstride1 = 1;
99      }          Nstride2 = N1;
100      else          mpi_info->split(NE0, &local_NE0, &e_offset0);
101      {          local_NE1 = NE1;
102          Nstride0 = N1 * N2;          e_offset1 = 0;
103          Nstride1 = 1;          local_NE2 = NE2;
104          Nstride2 = N1;          e_offset2 = 0;
105          mpi_info->split(NE0, &local_NE0, &e_offset0);      }
106          local_NE1 = NE1;      offset0 = e_offset0 * N_PER_E;
107          e_offset1 = 0;      offset1 = e_offset1 * N_PER_E;
108          local_NE2 = NE2;      offset2 = e_offset2 * N_PER_E;
109          e_offset2 = 0;      local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;
110      }      local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;
111      offset0 = e_offset0 * N_PER_E;      local_N2 = local_NE2 > 0 ? local_NE2 * N_PER_E + 1 : 0;
112      offset1 = e_offset1 * N_PER_E;  
113      offset2 = e_offset2 * N_PER_E;      /* get the number of surface elements */
114      local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;  
115      local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;      NFaceElements = 0;
116      local_N2 = local_NE2 > 0 ? local_NE2 * N_PER_E + 1 : 0;      if (local_NE2 > 0) {
117            NDOF2 = N2;
118      /* get the number of surface elements */          if (offset2 == 0)
119                NFaceElements += 2 * local_NE1 * local_NE0;     /* each face is split */
120      NFaceElements = 0;          if (local_NE2 + e_offset2 == NE2)
121      if (local_NE2 > 0)              NFaceElements += 2 * local_NE1 * local_NE0;
122      {      } else {
123          NDOF2 = N2;          NDOF2 = N2 - 1;
         if (offset2 == 0)  
         NFaceElements += 2 * local_NE1 * local_NE0; /* each face is split */  
         if (local_NE2 + e_offset2 == NE2)  
         NFaceElements += 2 * local_NE1 * local_NE0;  
     }  
     else  
     {  
         NDOF2 = N2 - 1;  
     }  
   
     if (local_NE0 > 0)  
     {  
         NDOF0 = N0;  
         if (e_offset0 == 0)  
         NFaceElements += 2 * local_NE1 * local_NE2;  
         if (local_NE0 + e_offset0 == NE0)  
         NFaceElements += 2 * local_NE1 * local_NE2;  
     }  
     else  
     {  
         NDOF0 = N0 - 1;  
     }  
   
     if (local_NE1 > 0)  
     {  
         NDOF1 = N1;  
         if (e_offset1 == 0)  
         NFaceElements += 2 * local_NE0 * local_NE2;  
         if (local_NE1 + e_offset1 == NE1)  
         NFaceElements += 2 * local_NE0 * local_NE2;  
     }  
     else  
     {  
         NDOF1 = N1 - 1;  
     }  
124      }      }
125    
126      /*  allocate tables: */      if (local_NE0 > 0) {
127      if (Dudley_noError())          NDOF0 = N0;
128      {          if (e_offset0 == 0)
129                NFaceElements += 2 * local_NE1 * local_NE2;
130            if (local_NE0 + e_offset0 == NE0)
131                NFaceElements += 2 * local_NE1 * local_NE2;
132        } else {
133            NDOF0 = N0 - 1;
134        }
135    
136      Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1 * local_N2);      if (local_NE1 > 0) {
137      /* we split the rectangular prism this code used to produce into 5 tetrahedrons */          NDOF1 = N1;
138      Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * local_NE2 * 5);          if (e_offset1 == 0)
139      /* each border face will be split in half */              NFaceElements += 2 * local_NE0 * local_NE2;
140      Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);          if (local_NE1 + e_offset1 == NE1)
141                NFaceElements += 2 * local_NE0 * local_NE2;
142        } else {
143            NDOF1 = N1 - 1;
144      }      }
145    
146      if (Dudley_noError())      /*  allocate tables: */
147      {      Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1 * local_N2);
148      int global_adjustment;      /* we split the rectangular prism this code used to produce into 5 tetrahedrons */
149        Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * local_NE2 * 5);
150        /* each border face will be split in half */
151        Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);
152    
153      /* create nodes */      int global_adjustment;
154    
155        /* create nodes */
156    
157  #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)  #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)
158      for (i2 = 0; i2 < local_N2; i2++)      for (i2 = 0; i2 < local_N2; i2++) {
159      {          for (i1 = 0; i1 < local_N1; i1++) {
160          for (i1 = 0; i1 < local_N1; i1++)              for (i0 = 0; i0 < local_N0; i0++) {
161          {                  k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2;
162          for (i0 = 0; i0 < local_N0; i0++)                  global_i0 = i0 + offset0;
163          {                  global_i1 = i1 + offset1;
164              k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2;                  global_i2 = i2 + offset2;
165              global_i0 = i0 + offset0;                  out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0];
166              global_i1 = i1 + offset1;                  out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];
167              global_i2 = i2 + offset2;                  out->Nodes->Coordinates[INDEX2(2, k, DIM)] = (double)global_i2 / (double)(N2 - 1) * Length[2];
168              out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0];                  out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2;
169              out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];                  out->Nodes->Tag[k] = 0;
170              out->Nodes->Coordinates[INDEX2(2, k, DIM)] = (double)global_i2 / (double)(N2 - 1) * Length[2];                  out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)
171              out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2;                      + Nstride1 * (global_i1 % NDOF1) + Nstride2 * (global_i2 % NDOF2);
172              out->Nodes->Tag[k] = 0;              }
173              out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)          }
174              + Nstride1 * (global_i1 % NDOF1) + Nstride2 * (global_i2 % NDOF2);      }
175          }      /*   set the elements: */
         }  
     }  
     /*   set the elements: */  
176    
177      global_adjustment = (offset0 + offset1 + offset2) % 2;  /* If we are not the only rank we may need to shift our pattern to match neighbours */      global_adjustment = (offset0 + offset1 + offset2) % 2;  /* If we are not the only rank we may need to shift our pattern to match neighbours */
178    
179      NN = out->Elements->numNodes;      NN = out->Elements->numNodes;
180  #pragma omp parallel for private(i0,i1,i2,k,node0)  #pragma omp parallel for private(i0,i1,i2,k,node0)
181      for (i2 = 0; i2 < local_NE2; i2++)      for (i2 = 0; i2 < local_NE2; i2++) {
182      {          for (i1 = 0; i1 < local_NE1; i1++) {
183          for (i1 = 0; i1 < local_NE1; i1++)              for (i0 = 0; i0 < local_NE0; i0++) {
184          {                  index_t res;
185          for (i0 = 0; i0 < local_NE0; i0++)                  index_t v[8];
186          {                  int j;
187              index_t res;                  k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2);
188              index_t v[8];                  node0 =
189              int j;                      Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
190              k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2);                      Nstride2 * N_PER_E * (i2 + e_offset2);
191              node0 =  
192              Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +                  res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2));
193              Nstride2 * N_PER_E * (i2 + e_offset2);                  for (j = 0; j < 5; ++j) {
194                        out->Elements->Id[k + j] = res + j;
195              res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2));                      out->Elements->Tag[k + j] = 0;
196              for (j = 0; j < 5; ++j)                      out->Elements->Owner[k + j] = myRank;
197              {                  }
198              out->Elements->Id[k + j] = res + j;  
199              out->Elements->Tag[k + j] = 0;  /*         in non-rotated orientation the points are numbered as follows:
200              out->Elements->Owner[k + j] = myRank;         The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/
201              }                  if ((global_adjustment + i0 + i1 + i2) % 2 == 0) {
202                        v[0] = node0;
203  /*     in non-rotated orientation the points are numbered as follows:                      v[1] = node0 + Nstride0;
204         The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/                      v[2] = node0 + Nstride1;
205              if ((global_adjustment + i0 + i1 + i2) % 2 == 0)                      v[3] = node0 + Nstride1 + Nstride0;
206              {                      v[4] = node0 + Nstride2;
207              v[0] = node0;                      v[5] = node0 + Nstride0 + Nstride2;
208              v[1] = node0 + Nstride0;                      v[6] = node0 + Nstride1 + Nstride2;
209              v[2] = node0 + Nstride1;                      v[7] = node0 + Nstride2 + Nstride1 + Nstride0;
210              v[3] = node0 + Nstride1 + Nstride0;                  } else {
211              v[4] = node0 + Nstride2;                      /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
212              v[5] = node0 + Nstride0 + Nstride2;  
213              v[6] = node0 + Nstride1 + Nstride2;                      v[0] = node0 + Nstride1;        /* node 0 ends up in position 2 */
214              v[7] = node0 + Nstride2 + Nstride1 + Nstride0;                      v[2] = node0 + Nstride1 + Nstride2;     /* node 2 ends up in position 6 */
215              }                      v[6] = node0 + Nstride2;        /* node 6 ends up in position 4 */
216              else                      v[4] = node0;   /* node 4 ends up in position 0 */
217              {  
218              /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                      v[1] = node0 + Nstride0 + Nstride1;     /* node 1 -> pos 3 */
219                        v[3] = node0 + Nstride2 + Nstride1 + Nstride0;  /* node 3-> pos 7 */
220              v[0] = node0 + Nstride1;    /* node 0 ends up in position 2 */                      v[7] = node0 + Nstride0 + Nstride2;     /* node 7 -> pos 5 */
221              v[2] = node0 + Nstride1 + Nstride2; /* node 2 ends up in position 6 */                      v[5] = node0 + Nstride0;        /* node 5 -> pos 1 */
222              v[6] = node0 + Nstride2;    /* node 6 ends up in position 4 */                  }
223              v[4] = node0;   /* node 4 ends up in position 0 */  
224                    /* elements nodes are numbered: centre, x, y, z */
225              v[1] = node0 + Nstride0 + Nstride1; /* node 1 -> pos 3 */  
226              v[3] = node0 + Nstride2 + Nstride1 + Nstride0;  /* node 3-> pos 7 */                  out->Elements->Nodes[INDEX2(0, k, NN)] = v[4];
227              v[7] = node0 + Nstride0 + Nstride2; /* node 7 -> pos 5 */                  out->Elements->Nodes[INDEX2(1, k, NN)] = v[5];
228              v[5] = node0 + Nstride0;    /* node 5 -> pos 1 */                  out->Elements->Nodes[INDEX2(2, k, NN)] = v[6];
229              }                  out->Elements->Nodes[INDEX2(3, k, NN)] = v[0];
230    
231              /* elements nodes are numbered: centre, x, y, z */                  out->Elements->Nodes[INDEX2(0, k + 1, NN)] = v[7];
232                    out->Elements->Nodes[INDEX2(1, k + 1, NN)] = v[6];
233              out->Elements->Nodes[INDEX2(0, k, NN)] = v[4];                  out->Elements->Nodes[INDEX2(2, k + 1, NN)] = v[5];
234              out->Elements->Nodes[INDEX2(1, k, NN)] = v[5];                  out->Elements->Nodes[INDEX2(3, k + 1, NN)] = v[3];
235              out->Elements->Nodes[INDEX2(2, k, NN)] = v[6];  
236              out->Elements->Nodes[INDEX2(3, k, NN)] = v[0];                  out->Elements->Nodes[INDEX2(0, k + 2, NN)] = v[2];
237                    out->Elements->Nodes[INDEX2(1, k + 2, NN)] = v[3];
238              out->Elements->Nodes[INDEX2(0, k + 1, NN)] = v[7];                  out->Elements->Nodes[INDEX2(2, k + 2, NN)] = v[0];
239              out->Elements->Nodes[INDEX2(1, k + 1, NN)] = v[6];                  out->Elements->Nodes[INDEX2(3, k + 2, NN)] = v[6];
240              out->Elements->Nodes[INDEX2(2, k + 1, NN)] = v[5];  
241              out->Elements->Nodes[INDEX2(3, k + 1, NN)] = v[3];                  out->Elements->Nodes[INDEX2(0, k + 3, NN)] = v[1];
242                    out->Elements->Nodes[INDEX2(1, k + 3, NN)] = v[0];
243              out->Elements->Nodes[INDEX2(0, k + 2, NN)] = v[2];                  out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3];
244              out->Elements->Nodes[INDEX2(1, k + 2, NN)] = v[3];                  out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5];
245              out->Elements->Nodes[INDEX2(2, k + 2, NN)] = v[0];  
246              out->Elements->Nodes[INDEX2(3, k + 2, NN)] = v[6];                  /* I can't work out where the center is for this one */
247                    out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5];
248              out->Elements->Nodes[INDEX2(0, k + 3, NN)] = v[1];                  out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0];
249              out->Elements->Nodes[INDEX2(1, k + 3, NN)] = v[0];                  out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6];
250              out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3];                  out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3];
251              out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5];              }
252            }
253              /* I can't work out where the center is for this one */      }                       /* end for */
254              out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5];      /* face elements */
255              out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0];      NN = out->FaceElements->numNodes;
256              out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6];      totalNECount = 5 * NE0 * NE1 * NE2;
257              out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3];      faceNECount = 0;
258          }      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
259        if (local_NE2 > 0) {
260          }          /* **  elements on boundary 100 (x3=0): */
261      }           /* end for */          if (e_offset2 == 0) {
     /* face elements */  
     NN = out->FaceElements->numNodes;  
     totalNECount = 5 * NE0 * NE1 * NE2;  
     faceNECount = 0;  
     /*   these are the quadrilateral elements on boundary 1 (x3=0): */  
     if (local_NE2 > 0)  
     {  
         /* **  elements on boundary 100 (x3=0): */  
         if (e_offset2 == 0)  
         {  
262  #pragma omp parallel for private(i0,i1,k,node0)  #pragma omp parallel for private(i0,i1,k,node0)
263          for (i1 = 0; i1 < local_NE1; i1++)              for (i1 = 0; i1 < local_NE1; i1++) {
264          {                  for (i0 = 0; i0 < local_NE0; i0++) {
265              for (i0 = 0; i0 < local_NE0; i0++)                      index_t res, n0, n1, n2, n3;
266              {                      k = 2 * (i0 + local_NE0 * i1) + faceNECount;
267              index_t res, n0, n1, n2, n3;                      node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);
268              k = 2 * (i0 + local_NE0 * i1) + faceNECount;                      res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
269              node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);                      out->FaceElements->Id[k] = res;
270              res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;                      out->FaceElements->Tag[k] = BOTTOMTAG;
271              out->FaceElements->Id[k] = res;                      out->FaceElements->Owner[k] = myRank;
272              out->FaceElements->Tag[k] = BOTTOMTAG;                      out->FaceElements->Id[k + 1] = res + 1;
273              out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Tag[k + 1] = BOTTOMTAG;
274              out->FaceElements->Id[k + 1] = res + 1;                      out->FaceElements->Owner[k + 1] = myRank;
275              out->FaceElements->Tag[k + 1] = BOTTOMTAG;  
276              out->FaceElements->Owner[k + 1] = myRank;                      n0 = node0;
277                        n1 = node0 + Nstride0;
278              n0 = node0;                      n2 = node0 + Nstride1;
279              n1 = node0 + Nstride0;                      n3 = node0 + Nstride0 + Nstride1;
280              n2 = node0 + Nstride1;  
281              n3 = node0 + Nstride0 + Nstride1;                      if ((global_adjustment + i0 + i1) % 2 == 0) {
282                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
283              if ((global_adjustment + i0 + i1) % 2 == 0)                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;
284              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;
285                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;  
286                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
287                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;
288                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
289                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;  
290                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;                      } else {
291                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
292                            out->FaceElements->Nodes[INDEX2(1, k, NN)] = n2;
293              }                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;
294              else  
295              {                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
296                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;
297                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
298                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n1;  
299                        }
300                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;                  }
301                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n2;              }
302                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;              faceNECount += 2 * local_NE1 * local_NE0;
303            }
304              }          totalNECount += 2 * NE1 * NE0;
305              }          /* **  elements on boundary 200 (x3=1) - Top */
306          }          if (local_NE2 + e_offset2 == NE2) {
         faceNECount += 2 * local_NE1 * local_NE0;  
         }  
         totalNECount += 2 * NE1 * NE0;  
         /* **  elements on boundary 200 (x3=1) - Top */  
         if (local_NE2 + e_offset2 == NE2)  
         {  
307  #pragma omp parallel for private(i0,i1,k,node0)  #pragma omp parallel for private(i0,i1,k,node0)
308          for (i1 = 0; i1 < local_NE1; i1++)              for (i1 = 0; i1 < local_NE1; i1++) {
309          {                  for (i0 = 0; i0 < local_NE0; i0++) {
             for (i0 = 0; i0 < local_NE0; i0++)  
             {  
   
             index_t res, n4, n5, n6, n7;  
             k = 2 * (i0 + local_NE0 * i1) + faceNECount;  
             node0 =  
                 Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +  
                 Nstride2 * N_PER_E * (NE2 - 1);  
   
             res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;  
             out->FaceElements->Id[k] = res;  
             out->FaceElements->Tag[k] = TOPTAG;  
             out->FaceElements->Owner[k] = myRank;  
             out->FaceElements->Id[k + 1] = res + 1;  
             out->FaceElements->Tag[k + 1] = TOPTAG;  
             out->FaceElements->Owner[k + 1] = myRank;  
   
             n4 = node0 + Nstride2;  
             n5 = node0 + Nstride0 + Nstride2;  
             n6 = node0 + Nstride1 + Nstride2;  
             n7 = node0 + Nstride1 + Nstride0 + Nstride2;  
   
             if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0)  
             {  
                 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;  
                 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;  
                 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;  
   
                 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n5;  
                 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;  
                 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;  
             }  
             else  
             {  
   
                 out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;  
                 out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;  
                 out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;  
   
                 out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;  
                 out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;  
                 out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;  
             }  
             }  
         }  
         faceNECount += 2 * local_NE1 * local_NE0;  
         }  
         totalNECount += 2 * NE1 * NE0;  
     }  
     if (local_NE0 > 0)  
     {  
         /* **  elements on boundary 001 (x1=0): - Left */  
310    
311          if (e_offset0 == 0)                      index_t res, n4, n5, n6, n7;
312          {                      k = 2 * (i0 + local_NE0 * i1) + faceNECount;
313                        node0 =
314                            Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
315                            Nstride2 * N_PER_E * (NE2 - 1);
316    
317                        res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;
318                        out->FaceElements->Id[k] = res;
319                        out->FaceElements->Tag[k] = TOPTAG;
320                        out->FaceElements->Owner[k] = myRank;
321                        out->FaceElements->Id[k + 1] = res + 1;
322                        out->FaceElements->Tag[k + 1] = TOPTAG;
323                        out->FaceElements->Owner[k + 1] = myRank;
324    
325                        n4 = node0 + Nstride2;
326                        n5 = node0 + Nstride0 + Nstride2;
327                        n6 = node0 + Nstride1 + Nstride2;
328                        n7 = node0 + Nstride1 + Nstride0 + Nstride2;
329    
330                        if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0) {
331                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;
332                            out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;
333                            out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;
334    
335                            out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n5;
336                            out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
337                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;
338                        } else {
339    
340                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;
341                            out->FaceElements->Nodes[INDEX2(1, k, NN)] = n5;
342                            out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;
343    
344                            out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;
345                            out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
346                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n6;
347                        }
348                    }
349                }
350                faceNECount += 2 * local_NE1 * local_NE0;
351            }
352            totalNECount += 2 * NE1 * NE0;
353        }
354        if (local_NE0 > 0) {
355            /* **  elements on boundary 001 (x1=0): - Left */
356    
357            if (e_offset0 == 0) {
358  #pragma omp parallel for private(i1,i2,k,node0)  #pragma omp parallel for private(i1,i2,k,node0)
359          for (i2 = 0; i2 < local_NE2; i2++)              for (i2 = 0; i2 < local_NE2; i2++) {
360          {                  for (i1 = 0; i1 < local_NE1; i1++) {
361              for (i1 = 0; i1 < local_NE1; i1++)  
362              {                      index_t res, n0, n2, n4, n6;
363                        k = 2 * (i1 + local_NE1 * i2) + faceNECount;
364              index_t res, n0, n2, n4, n6;                      node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2);
365              k = 2 * (i1 + local_NE1 * i2) + faceNECount;                      res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
366              node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2);                      out->FaceElements->Id[k] = res;
367              res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;                      out->FaceElements->Tag[k] = LEFTTAG;
368              out->FaceElements->Id[k] = res;                      out->FaceElements->Owner[k] = myRank;
369              out->FaceElements->Tag[k] = LEFTTAG;                      out->FaceElements->Id[k + 1] = res + 1;
370              out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Tag[k + 1] = LEFTTAG;
371              out->FaceElements->Id[k + 1] = res + 1;                      out->FaceElements->Owner[k + 1] = myRank;
372              out->FaceElements->Tag[k + 1] = LEFTTAG;  
373              out->FaceElements->Owner[k + 1] = myRank;                      n0 = node0;
374                        n2 = node0 + Nstride1;
375              n0 = node0;                      n4 = node0 + Nstride2;
376              n2 = node0 + Nstride1;                      n6 = node0 + Nstride1 + Nstride2;
377              n4 = node0 + Nstride2;  
378              n6 = node0 + Nstride1 + Nstride2;                      if ((global_adjustment + 0 + i1 + i2) % 2 == 0) {
379                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
380              if ((global_adjustment + 0 + i1 + i2) % 2 == 0)                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;
381              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;
382                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;  
383                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
384                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n6;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;
385                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;
386                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;                      } else {
387                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
388                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;  
389              }                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
390              else                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;
391              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2;
392                  /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */  
393                            out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;
394                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;
395                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;
396                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2;                      }
397                    }
398                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n4;              }
399                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;              faceNECount += 2 * local_NE1 * local_NE2;
400                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;          }
401              }          totalNECount += 2 * NE1 * NE2;
402              }          /* **  elements on boundary 002 (x1=1): - Right */
403          }          if (local_NE0 + e_offset0 == NE0) {
         faceNECount += 2 * local_NE1 * local_NE2;  
         }  
         totalNECount += 2 * NE1 * NE2;  
         /* **  elements on boundary 002 (x1=1): - Right */  
         if (local_NE0 + e_offset0 == NE0)  
         {  
404  #pragma omp parallel for private(i1,i2,k,node0)  #pragma omp parallel for private(i1,i2,k,node0)
405          for (i2 = 0; i2 < local_NE2; i2++)              for (i2 = 0; i2 < local_NE2; i2++) {
406          {                  for (i1 = 0; i1 < local_NE1; i1++) {
407              for (i1 = 0; i1 < local_NE1; i1++)                      index_t res, n1, n3, n5, n7;
408              {                      k = 2 * (i1 + local_NE1 * i2) + faceNECount;
409              index_t res, n1, n3, n5, n7;  
410              k = 2 * (i1 + local_NE1 * i2) + faceNECount;                      node0 =
411                            Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) +
412              node0 =                          Nstride2 * N_PER_E * (i2 + e_offset2);
413                  Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) +                      res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
414                  Nstride2 * N_PER_E * (i2 + e_offset2);                      out->FaceElements->Id[k] = res;
415              res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;                      out->FaceElements->Tag[k] = RIGHTTAG;
416              out->FaceElements->Id[k] = res;                      out->FaceElements->Owner[k] = myRank;
417              out->FaceElements->Tag[k] = RIGHTTAG;                      out->FaceElements->Id[k + 1] = res + 1;
418              out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Tag[k + 1] = RIGHTTAG;
419              out->FaceElements->Id[k + 1] = res + 1;                      out->FaceElements->Owner[k + 1] = myRank;
420              out->FaceElements->Tag[k + 1] = RIGHTTAG;  
421              out->FaceElements->Owner[k + 1] = myRank;                      n1 = node0 + Nstride0;
422                        n3 = node0 + Nstride0 + Nstride1;
423              n1 = node0 + Nstride0;                      n5 = node0 + Nstride0 + Nstride2;
424              n3 = node0 + Nstride0 + Nstride1;                      n7 = node0 + Nstride0 + Nstride1 + Nstride2;
425              n5 = node0 + Nstride0 + Nstride2;                      if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0) {
426              n7 = node0 + Nstride0 + Nstride1 + Nstride2;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
427              if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0)                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;
428              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
429                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;  
430                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n3;
431                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
432                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5;
433                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n3;                      } else {
434                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
435                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5;  
436              }                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
437              else                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7;
438              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
439                  /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */  
440                            out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
441                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n3;
442                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n7;
443                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;                      }
444                    }
445                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;              }
446                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n3;              faceNECount += 2 * local_NE1 * local_NE2;
447                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n7;          }
448              }          totalNECount += 2 * NE1 * NE2;
449              }      }
450          }      if (local_NE1 > 0) {
451          faceNECount += 2 * local_NE1 * local_NE2;          /* **  elements on boundary 010 (x2=0): -Front */
452          }          if (e_offset1 == 0) {
         totalNECount += 2 * NE1 * NE2;  
     }  
     if (local_NE1 > 0)  
     {  
         /* **  elements on boundary 010 (x2=0): -Front */  
         if (e_offset1 == 0)  
         {  
453  #pragma omp parallel for private(i0,i2,k,node0)  #pragma omp parallel for private(i0,i2,k,node0)
454          for (i2 = 0; i2 < local_NE2; i2++)              for (i2 = 0; i2 < local_NE2; i2++) {
455          {                  for (i0 = 0; i0 < local_NE0; i0++) {
456              for (i0 = 0; i0 < local_NE0; i0++)                      index_t res, n0, n1, n4, n5;
457              {                      k = 2 * (i0 + local_NE0 * i2) + faceNECount;
458              index_t res, n0, n1, n4, n5;                      node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2);
459              k = 2 * (i0 + local_NE0 * i2) + faceNECount;                      res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount;
460              node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2);                      out->FaceElements->Id[k] = res;
461              res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount;                      out->FaceElements->Tag[k] = FRONTTAG;
462              out->FaceElements->Id[k] = res;                      out->FaceElements->Owner[k] = myRank;
463              out->FaceElements->Tag[k] = FRONTTAG;                      out->FaceElements->Id[k + 1] = res + 1;
464              out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Tag[k + 1] = FRONTTAG;
465              out->FaceElements->Id[k + 1] = res + 1;                      out->FaceElements->Owner[k + 1] = myRank;
466              out->FaceElements->Tag[k + 1] = FRONTTAG;  
467              out->FaceElements->Owner[k + 1] = myRank;                      n0 = node0;
468                        n1 = node0 + Nstride0;
469              n0 = node0;                      n4 = node0 + Nstride2;
470              n1 = node0 + Nstride0;                      n5 = node0 + Nstride0 + Nstride2;
471              n4 = node0 + Nstride2;  
472              n5 = node0 + Nstride0 + Nstride2;                      if ((global_adjustment + i0 + 0 + i2) % 2 == 0) {
473                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
474              if ((global_adjustment + i0 + 0 + i2) % 2 == 0)                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;
475              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
476                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;  
477                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;
478                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;
479                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;
480                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n0;  
481                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;                      } else {
482                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
483    
484              }                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
485              else                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;
486              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4;
487                  /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */  
488                            out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;
489                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;
490                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;
491                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4;  
492                        }
493                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n1;                  }
494                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n5;              }
495                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;              faceNECount += 2 * local_NE0 * local_NE2;
496            }
497              }          totalNECount += 2 * NE0 * NE2;
498              }          /* **  elements on boundary 020 (x2=1): - Back */
499          }          if (local_NE1 + e_offset1 == NE1) {
         faceNECount += 2 * local_NE0 * local_NE2;  
         }  
         totalNECount += 2 * NE0 * NE2;  
         /* **  elements on boundary 020 (x2=1): - Back */  
         if (local_NE1 + e_offset1 == NE1)  
         {  
500  #pragma omp parallel for private(i0,i2,k,node0)  #pragma omp parallel for private(i0,i2,k,node0)
501          for (i2 = 0; i2 < local_NE2; i2++)              for (i2 = 0; i2 < local_NE2; i2++) {
502          {                  for (i0 = 0; i0 < local_NE0; i0++) {
503              for (i0 = 0; i0 < local_NE0; i0++)                      index_t res, n2, n6, n7, n3;
504              {                      k = 2 * (i0 + local_NE0 * i2) + faceNECount;
505              index_t res, n2, n6, n7, n3;                      node0 =
506              k = 2 * (i0 + local_NE0 * i2) + faceNECount;                          Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1) +
507              node0 =                          Nstride2 * N_PER_E * (i2 + e_offset2);
508                  Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1) +                      res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount;
509                  Nstride2 * N_PER_E * (i2 + e_offset2);                      out->FaceElements->Id[k] = res;
510              res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount;                      out->FaceElements->Tag[k] = BACKTAG;
511              out->FaceElements->Id[k] = res;                      out->FaceElements->Owner[k] = myRank;
512              out->FaceElements->Tag[k] = BACKTAG;                      out->FaceElements->Id[k + 1] = res + 1;
513              out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Tag[k + 1] = BACKTAG;
514              out->FaceElements->Id[k + 1] = res + 1;                      out->FaceElements->Owner[k + 1] = myRank;
515              out->FaceElements->Tag[k + 1] = BACKTAG;  
516              out->FaceElements->Owner[k + 1] = myRank;                      n2 = node0 + Nstride1;
517                        n6 = node0 + Nstride1 + Nstride2;
518              n2 = node0 + Nstride1;                      n7 = node0 + Nstride0 + Nstride1 + Nstride2;
519              n6 = node0 + Nstride1 + Nstride2;                      n3 = node0 + Nstride0 + Nstride1;
520              n7 = node0 + Nstride0 + Nstride1 + Nstride2;  
521              n3 = node0 + Nstride0 + Nstride1;                      if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0) {
522                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
523              if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0)                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;
524              {                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n3;
525                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;  
526                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n6;
527                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
528                            out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
529                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n6;  
530                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;                      } else {
531                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */
532                            out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
533              }                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;
534              else                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;
535              {  
536                  /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2;
537                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
538                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
539                  out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;  
540                        }
541                  out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2;                  }
542                  out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;              }
543                  out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;              faceNECount += 2 * local_NE0 * local_NE2;
544            }
545              }          totalNECount += 2 * NE0 * NE2;
             }  
         }  
         faceNECount += 2 * local_NE0 * local_NE2;  
         }  
         totalNECount += 2 * NE0 * NE2;  
     }  
     }  
     if (Dudley_noError())  
     {  
     /* add tag names */  
     Dudley_Mesh_addTagMap(out, "top", TOPTAG);  
     Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG);  
     Dudley_Mesh_addTagMap(out, "left", LEFTTAG);  
     Dudley_Mesh_addTagMap(out, "right", RIGHTTAG);  
     Dudley_Mesh_addTagMap(out, "front", FRONTTAG);  
     Dudley_Mesh_addTagMap(out, "back", BACKTAG);  
     }  
     /* prepare mesh for further calculations: */  
     if (Dudley_noError())  
     {  
     Dudley_Mesh_resolveNodeIds(out);  
     }  
     if (Dudley_noError())  
     {  
     Dudley_Mesh_prepare(out, optimize);  
546      }      }
547        /* add tag names */
548        Dudley_Mesh_addTagMap(out, "top", TOPTAG);
549        Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG);
550        Dudley_Mesh_addTagMap(out, "left", LEFTTAG);
551        Dudley_Mesh_addTagMap(out, "right", RIGHTTAG);
552        Dudley_Mesh_addTagMap(out, "front", FRONTTAG);
553        Dudley_Mesh_addTagMap(out, "back", BACKTAG);
554        // prepare mesh for further calculations
555        Dudley_Mesh_resolveNodeIds(out);
556        Dudley_Mesh_prepare(out, optimize);
557    
     if (!Dudley_noError())  
     {  
     Dudley_Mesh_free(out);  
     }  
     /* free up memory */  
558      return out;      return out;
559  }  }
560    
561    } // namespace dudley
562    

Legend:
Removed from v.6008  
changed lines
  Added in v.6009

  ViewVC Help
Powered by ViewVC 1.1.26