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

revision 6078 by caltinay, Wed Mar 2 04:13:26 2016 UTC revision 6079 by caltinay, Mon Mar 21 12:22:38 2016 UTC
# Line 14  Line 14
14  *  *
15  *****************************************************************************/  *****************************************************************************/
16
/****************************************************************************/

/*   Dudley: generates rectangular meshes */

/*   Generates a numElements[0] x numElements[1] x numElements[2] mesh with first order elements (Hex8) in the brick */
/*   [0,Length[0]] x [0,Length[1]] x [0,Length[2]]. order is the desired accuracy of the */
/*   integration scheme. */

/****************************************************************************/

17  #include "TriangularMesh.h"  #include "TriangularMesh.h"
18
19    using escript::DataTypes::real_t;
20
21  #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_))  #define MAX3(_arg1_,_arg2_,_arg3_) std::max(_arg1_,std::max(_arg2_,_arg3_))
22
23  namespace dudley {  namespace dudley {
24
25  // Be careful reading this function. The X? and NStride? are 1,2,3  // Be careful reading this function. The X? and NStride? are 1,2,3
26  // but the loop vars are 0,1,2  // but the loop vars are 0,1,2
27  Dudley_Mesh *Dudley_TriangularMesh_Tet4(dim_t * numElements,  Mesh* TriangularMesh_Tet4(const dim_t* numElements, const double* length,
28                                          double *Length, index_t order, index_t reduced_order, bool optimize, escript::JMPI& mpi_info)                            bool optimize, escript::JMPI mpiInfo)
29  {  {
30  #define N_PER_E 1      const int DIM = 3;
#define DIM 3
dim_t N0, N1, N2, NE0, NE1, NE2, i0, i1, i2, k, Nstride0 = 0, Nstride1 = 0, Nstride2 =
0, local_NE0, local_NE1, local_NE2, local_N0 = 0, local_N1 = 0, local_N2 = 0;
dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NDOF2 = 0, NFaceElements = 0, NN;
index_t node0, myRank, e_offset2, e_offset1, e_offset0 = 0, offset1 = 0, offset2 = 0, offset0 =
0, global_i0, global_i1, global_i2;
Dudley_Mesh *out;
char name[50];
31  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
32      double time0 = Dudley_timer();      double time0 = Dudley_timer();
33  #endif  #endif
# Line 55  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 39  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
39      const int FRONTTAG = 10;    /* boundary x2=0 */      const int FRONTTAG = 10;    /* boundary x2=0 */
40      const int BACKTAG = 20;     /* boundary x2=1 */      const int BACKTAG = 20;     /* boundary x2=1 */
41
42      /* get MPI information */      const int myRank = mpiInfo->rank;
myRank = mpi_info->rank;

/* set up the global dimensions of the mesh */
NE0 = std::max(1, numElements[0]);
NE1 = std::max(1, numElements[1]);
NE2 = std::max(1, numElements[2]);
N0 = N_PER_E * NE0 + 1;
N1 = N_PER_E * NE1 + 1;
N2 = N_PER_E * NE2 + 1;

/*  allocate mesh: */
sprintf(name, "Triangular %d x %d x %d (x 5) mesh", N0, N1, N2);
out = Dudley_Mesh_alloc(name, DIM, mpi_info);

Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));
Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));
Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tet4, mpi_info));
43
44      /* work out the largest dimension */      // set up the global dimensions of the mesh
45        const dim_t NE0 = std::max(dim_t(1), numElements[0]);
46        const dim_t NE1 = std::max(dim_t(1), numElements[1]);
47        const dim_t NE2 = std::max(dim_t(1), numElements[2]);
48        const dim_t N0 = NE0 + 1;
49        const dim_t N1 = NE1 + 1;
50        const dim_t N2 = NE2 + 1;
51
52        // allocate mesh
53        std::stringstream name;
54        name << "Triangular " << N0 << " x " << N1 << " x " << N2 << " (x 5) mesh";
55        Mesh* out = new Mesh(name.str(), DIM, mpiInfo);
56
57        out->setPoints(new ElementFile(Dudley_Point1, mpiInfo));
58        out->setFaceElements(new ElementFile(Dudley_Tri3, mpiInfo));
59        out->setElements(new ElementFile(Dudley_Tet4, mpiInfo));
60
61        dim_t Nstride0, Nstride1, Nstride2;
62        dim_t local_NE0, local_NE1, local_NE2;
63        index_t e_offset0, e_offset1, e_offset2;
64        // work out the largest dimension
65      if (N2 == MAX3(N0, N1, N2)) {      if (N2 == MAX3(N0, N1, N2)) {
66          Nstride0 = 1;          Nstride0 = 1;
67          Nstride1 = N0;          Nstride1 = N0;
# Line 83  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 70  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
70          e_offset0 = 0;          e_offset0 = 0;
71          local_NE1 = NE1;          local_NE1 = NE1;
72          e_offset1 = 0;          e_offset1 = 0;
73          mpi_info->split(NE2, &local_NE2, &e_offset2);          mpiInfo->split(NE2, &local_NE2, &e_offset2);
74      } else if (N1 == MAX3(N0, N1, N2)) {      } else if (N1 == MAX3(N0, N1, N2)) {
75          Nstride0 = N2;          Nstride0 = N2;
76          Nstride1 = N0 * N2;          Nstride1 = N0 * N2;
77          Nstride2 = 1;          Nstride2 = 1;
78          local_NE0 = NE0;          local_NE0 = NE0;
79          e_offset0 = 0;          e_offset0 = 0;
80          mpi_info->split(NE1, &local_NE1, &e_offset1);          mpiInfo->split(NE1, &local_NE1, &e_offset1);
81          local_NE2 = NE2;          local_NE2 = NE2;
82          e_offset2 = 0;          e_offset2 = 0;
83      } else {      } else {
84          Nstride0 = N1 * N2;          Nstride0 = N1 * N2;
85          Nstride1 = 1;          Nstride1 = 1;
86          Nstride2 = N1;          Nstride2 = N1;
87          mpi_info->split(NE0, &local_NE0, &e_offset0);          mpiInfo->split(NE0, &local_NE0, &e_offset0);
88          local_NE1 = NE1;          local_NE1 = NE1;
89          e_offset1 = 0;          e_offset1 = 0;
90          local_NE2 = NE2;          local_NE2 = NE2;
91          e_offset2 = 0;          e_offset2 = 0;
92      }      }
93      offset0 = e_offset0 * N_PER_E;      const index_t offset0 = e_offset0;
94      offset1 = e_offset1 * N_PER_E;      const index_t offset1 = e_offset1;
95      offset2 = e_offset2 * N_PER_E;      const index_t offset2 = e_offset2;
96      local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;      const dim_t local_N0 = local_NE0 > 0 ? local_NE0 + 1 : 0;
97      local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;      const dim_t local_N1 = local_NE1 > 0 ? local_NE1 + 1 : 0;
98      local_N2 = local_NE2 > 0 ? local_NE2 * N_PER_E + 1 : 0;      const dim_t local_N2 = local_NE2 > 0 ? local_NE2 + 1 : 0;
99
100      /* get the number of surface elements */      // get the number of surface elements
101        dim_t NFaceElements = 0;
102      NFaceElements = 0;      dim_t NDOF0, NDOF1, NDOF2;
103      if (local_NE2 > 0) {      if (local_NE2 > 0) {
104          NDOF2 = N2;          NDOF2 = N2;
105          if (offset2 == 0)          if (offset2 == 0)
106              NFaceElements += 2 * local_NE1 * local_NE0;     /* each face is split */              NFaceElements += 2 * local_NE1 * local_NE0; // each face is split
107          if (local_NE2 + e_offset2 == NE2)          if (local_NE2 + e_offset2 == NE2)
108              NFaceElements += 2 * local_NE1 * local_NE0;              NFaceElements += 2 * local_NE1 * local_NE0;
109      } else {      } else {
# Line 143  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 130  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
130          NDOF1 = N1 - 1;          NDOF1 = N1 - 1;
131      }      }
132
133      /*  allocate tables: */      // allocate tables
134      Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1 * local_N2);      out->Nodes->allocTable(local_N0 * local_N1 * local_N2);
135      /* we split the rectangular prism this code used to produce into 5 tetrahedrons */      // we split the rectangular prism this code used to produce into 5
136      Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * local_NE2 * 5);      // tetrahedra
137      /* each border face will be split in half */      out->Elements->allocTable(local_NE0 * local_NE1 * local_NE2 * 5);
138      Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);      // each border face will be split in half
139        out->FaceElements->allocTable(NFaceElements);
141        // create nodes
142      /* create nodes */  #pragma omp parallel for
143        for (index_t i2 = 0; i2 < local_N2; i2++) {
144  #pragma omp parallel for private(i0,i1,i2,k,global_i0,global_i1,global_i2)          for (index_t i1 = 0; i1 < local_N1; i1++) {
145      for (i2 = 0; i2 < local_N2; i2++) {              for (index_t i0 = 0; i0 < local_N0; i0++) {
146          for (i1 = 0; i1 < local_N1; i1++) {                  const index_t k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2;
147              for (i0 = 0; i0 < local_N0; i0++) {                  const index_t global_i0 = i0 + offset0;
148                  k = i0 + local_N0 * i1 + local_N0 * local_N1 * i2;                  const index_t global_i1 = i1 + offset1;
149                  global_i0 = i0 + offset0;                  const index_t global_i2 = i2 + offset2;
150                  global_i1 = i1 + offset1;                  out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (real_t)global_i0 / (real_t)(N0 - 1) * length[0];
151                  global_i2 = i2 + offset2;                  out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (real_t)global_i1 / (real_t)(N1 - 1) * length[1];
152                  out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0];                  out->Nodes->Coordinates[INDEX2(2, k, DIM)] = (real_t)global_i2 / (real_t)(N2 - 1) * length[2];
out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];
out->Nodes->Coordinates[INDEX2(2, k, DIM)] = (double)global_i2 / (double)(N2 - 1) * Length[2];
153                  out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2;                  out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1 + Nstride2 * global_i2;
154                  out->Nodes->Tag[k] = 0;                  out->Nodes->Tag[k] = 0;
155                  out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0)                  out->Nodes->globalDegreesOfFreedom[k] =
156                      + Nstride1 * (global_i1 % NDOF1) + Nstride2 * (global_i2 % NDOF2);                                      Nstride0 * (global_i0 % NDOF0)
157                                        + Nstride1 * (global_i1 % NDOF1)
158                                        + Nstride2 * (global_i2 % NDOF2);
159              }              }
160          }          }
161      }      }
/*   set the elements: */

global_adjustment = (offset0 + offset1 + offset2) % 2;  /* If we are not the only rank we may need to shift our pattern to match neighbours */

NN = out->Elements->numNodes;
#pragma omp parallel for private(i0,i1,i2,k,node0)
for (i2 = 0; i2 < local_NE2; i2++) {
for (i1 = 0; i1 < local_NE1; i1++) {
for (i0 = 0; i0 < local_NE0; i0++) {
index_t res;
index_t v[8];
int j;
k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2);
node0 =
Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
Nstride2 * N_PER_E * (i2 + e_offset2);
162
163                  res = 5 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1) + NE0 * NE1 * (i2 + e_offset2));      // set the elements
164                  for (j = 0; j < 5; ++j) {      // If we are not the only rank we may need to shift our pattern to match
165        // neighbours
166        int global_adjustment = (offset0 + offset1 + offset2) % 2;
167
168        int NN = out->Elements->numNodes;
169    #pragma omp parallel for
170        for (index_t i2 = 0; i2 < local_NE2; i2++) {
171            for (index_t i1 = 0; i1 < local_NE1; i1++) {
172                for (index_t i0 = 0; i0 < local_NE0; i0++) {
173                    const index_t k = 5 * (i0 + local_NE0 * i1 + local_NE0 * local_NE1 * i2);
174                    const index_t node0 = Nstride0 * (i0 + e_offset0)
175                                        + Nstride1 * (i1 + e_offset1)
176                                        + Nstride2 * (i2 + e_offset2);
177
178                    const index_t res = 5 * ((i0 + e_offset0)
179                                      + NE0 * (i1 + e_offset1)
180                                      + NE0 * NE1 * (i2 + e_offset2));
181                    for (int j = 0; j < 5; ++j) {
182                      out->Elements->Id[k + j] = res + j;                      out->Elements->Id[k + j] = res + j;
183                      out->Elements->Tag[k + j] = 0;                      out->Elements->Tag[k + j] = 0;
184                      out->Elements->Owner[k + j] = myRank;                      out->Elements->Owner[k + j] = myRank;
185                  }                  }
186
187  /*         in non-rotated orientation the points are numbered as follows:                  // in non-rotated orientation the points are numbered as
188         The bottom face (anticlockwise= 0,1,3,2), top face (anticlockwise 4,5,7,6)*/                  // follows:
189                    // The bottom face (anticlockwise = 0,1,3,2),
190                    // top face (anticlockwise 4,5,7,6)
191                    index_t v[8];
192                  if ((global_adjustment + i0 + i1 + i2) % 2 == 0) {                  if ((global_adjustment + i0 + i1 + i2) % 2 == 0) {
193                      v[0] = node0;                      v[0] = node0;
194                      v[1] = node0 + Nstride0;                      v[1] = node0 + Nstride0;
# Line 208  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 199  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
199                      v[6] = node0 + Nstride1 + Nstride2;                      v[6] = node0 + Nstride1 + Nstride2;
200                      v[7] = node0 + Nstride2 + Nstride1 + Nstride0;                      v[7] = node0 + Nstride2 + Nstride1 + Nstride0;
201                  } else {                  } else {
202                      /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                      // this form is rotated around the 0,2,4,6 face clockwise
203                        // 90 degrees
204                      v[0] = node0 + Nstride1;        /* node 0 ends up in position 2 */                      v[0] = node0 + Nstride1; // node 0 ends up in position 2
205                      v[2] = node0 + Nstride1 + Nstride2;     /* node 2 ends up in position 6 */                      v[2] = node0 + Nstride1 + Nstride2; // node 2 ends up in position 6
206                      v[6] = node0 + Nstride2;        /* node 6 ends up in position 4 */                      v[6] = node0 + Nstride2; // node 6 ends up in position 4
207                      v[4] = node0;   /* node 4 ends up in position 0 */                      v[4] = node0; // node 4 ends up in position 0
208                        v[1] = node0 + Nstride0 + Nstride1; // node 1 -> pos 3
209                      v[1] = node0 + Nstride0 + Nstride1;     /* node 1 -> pos 3 */                      v[3] = node0 + Nstride2 + Nstride1 + Nstride0; // node 3 -> pos 7
210                      v[3] = node0 + Nstride2 + Nstride1 + Nstride0;  /* node 3-> pos 7 */                      v[7] = node0 + Nstride0 + Nstride2; // node 7 -> pos 5
211                      v[7] = node0 + Nstride0 + Nstride2;     /* node 7 -> pos 5 */                      v[5] = node0 + Nstride0; // node 5 -> pos 1
v[5] = node0 + Nstride0;        /* node 5 -> pos 1 */
212                  }                  }
213
214                  /* elements nodes are numbered: centre, x, y, z */                  // elements nodes are numbered: centre, x, y, z

215                  out->Elements->Nodes[INDEX2(0, k, NN)] = v[4];                  out->Elements->Nodes[INDEX2(0, k, NN)] = v[4];
216                  out->Elements->Nodes[INDEX2(1, k, NN)] = v[5];                  out->Elements->Nodes[INDEX2(1, k, NN)] = v[5];
217                  out->Elements->Nodes[INDEX2(2, k, NN)] = v[6];                  out->Elements->Nodes[INDEX2(2, k, NN)] = v[6];
# Line 243  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 232  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
232                  out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3];                  out->Elements->Nodes[INDEX2(2, k + 3, NN)] = v[3];
233                  out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5];                  out->Elements->Nodes[INDEX2(3, k + 3, NN)] = v[5];
234
235                  /* I can't work out where the center is for this one */                  // I can't work out where the center is for this one
236                  out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5];                  out->Elements->Nodes[INDEX2(0, k + 4, NN)] = v[5];
237                  out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0];                  out->Elements->Nodes[INDEX2(1, k + 4, NN)] = v[0];
238                  out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6];                  out->Elements->Nodes[INDEX2(2, k + 4, NN)] = v[6];
239                  out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3];                  out->Elements->Nodes[INDEX2(3, k + 4, NN)] = v[3];
240              }              }
241          }          }
242      }                       /* end for */      } // for all elements
243      /* face elements */
244        // face elements
245      NN = out->FaceElements->numNodes;      NN = out->FaceElements->numNodes;
246      totalNECount = 5 * NE0 * NE1 * NE2;      dim_t totalNECount = 5 * NE0 * NE1 * NE2;
247      faceNECount = 0;      dim_t faceNECount = 0;
248      /*   these are the quadrilateral elements on boundary 1 (x3=0): */
249        // these are the quadrilateral elements on boundary 1 (x3=0)
250      if (local_NE2 > 0) {      if (local_NE2 > 0) {
251          /* **  elements on boundary 100 (x3=0): */          // ** elements on boundary 100 (x3=0)
252          if (e_offset2 == 0) {          if (e_offset2 == 0) {
253  #pragma omp parallel for private(i0,i1,k,node0)  #pragma omp parallel for
254              for (i1 = 0; i1 < local_NE1; i1++) {              for (index_t i1 = 0; i1 < local_NE1; i1++) {
255                  for (i0 = 0; i0 < local_NE0; i0++) {                  for (index_t i0 = 0; i0 < local_NE0; i0++) {
256                      index_t res, n0, n1, n2, n3;                      const index_t k = 2 * (i0 + local_NE0 * i1) + faceNECount;
257                      k = 2 * (i0 + local_NE0 * i1) + faceNECount;                      const index_t node0 = Nstride0 * (i0 + e_offset0)
258                      node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);                                          + Nstride1 * (i1 + e_offset1);
259                      res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;                      const index_t res = 2 * (i0 + e_offset0)
260                                         + NE0 * (i1 + e_offset1) + totalNECount;
261                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
262                      out->FaceElements->Tag[k] = BOTTOMTAG;                      out->FaceElements->Tag[k] = BOTTOMTAG;
263                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 273  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 265  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
265                      out->FaceElements->Tag[k + 1] = BOTTOMTAG;                      out->FaceElements->Tag[k + 1] = BOTTOMTAG;
266                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
267
268                      n0 = node0;                      const index_t n0 = node0;
269                      n1 = node0 + Nstride0;                      const index_t n1 = node0 + Nstride0;
270                      n2 = node0 + Nstride1;                      const index_t n2 = node0 + Nstride1;
271                      n3 = node0 + Nstride0 + Nstride1;                      const index_t n3 = node0 + Nstride0 + Nstride1;
272
273                      if ((global_adjustment + i0 + i1) % 2 == 0) {                      if ((global_adjustment + i0 + i1) % 2 == 0) {
274                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
# Line 302  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 294  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
294              faceNECount += 2 * local_NE1 * local_NE0;              faceNECount += 2 * local_NE1 * local_NE0;
295          }          }
296          totalNECount += 2 * NE1 * NE0;          totalNECount += 2 * NE1 * NE0;
297          /* **  elements on boundary 200 (x3=1) - Top */          // ** elements on boundary 200 (x3=1) - Top
298          if (local_NE2 + e_offset2 == NE2) {          if (local_NE2 + e_offset2 == NE2) {
299  #pragma omp parallel for private(i0,i1,k,node0)  #pragma omp parallel for
300              for (i1 = 0; i1 < local_NE1; i1++) {              for (index_t i1 = 0; i1 < local_NE1; i1++) {
301                  for (i0 = 0; i0 < local_NE0; i0++) {                  for (index_t i0 = 0; i0 < local_NE0; i0++) {
302                        const index_t k = 2 * (i0 + local_NE0 * i1) + faceNECount;
303                      index_t res, n4, n5, n6, n7;                      const index_t node0 = Nstride0 * (i0 + e_offset0)
304                      k = 2 * (i0 + local_NE0 * i1) + faceNECount;                                          + Nstride1 * (i1 + e_offset1)
305                      node0 =                                          + Nstride2 * (NE2 - 1);
Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1) +
Nstride2 * N_PER_E * (NE2 - 1);
306
307                      res = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1)) + totalNECount;                      const index_t res = 2 * (i0 + e_offset0)
308                                          + NE0 * (i1 + e_offset1) + totalNECount;
309                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
310                      out->FaceElements->Tag[k] = TOPTAG;                      out->FaceElements->Tag[k] = TOPTAG;
311                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 322  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 313  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
313                      out->FaceElements->Tag[k + 1] = TOPTAG;                      out->FaceElements->Tag[k + 1] = TOPTAG;
314                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
315
316                      n4 = node0 + Nstride2;                      const index_t n4 = node0 + Nstride2;
317                      n5 = node0 + Nstride0 + Nstride2;                      const index_t n5 = node0 + Nstride0 + Nstride2;
318                      n6 = node0 + Nstride1 + Nstride2;                      const index_t n6 = node0 + Nstride1 + Nstride2;
319                      n7 = node0 + Nstride1 + Nstride0 + Nstride2;                      const index_t n7 = node0 + Nstride1 + Nstride0 + Nstride2;
320
321                      if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0) {                      if ((global_adjustment + i0 + i1 + local_NE2 - 1) % 2 == 0) {
322                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n4;
# Line 351  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 342  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
342          }          }
343          totalNECount += 2 * NE1 * NE0;          totalNECount += 2 * NE1 * NE0;
344      }      }
if (local_NE0 > 0) {
/* **  elements on boundary 001 (x1=0): - Left */
345
346        if (local_NE0 > 0) {
347            // ** elements on boundary 001 (x1=0) - Left
348          if (e_offset0 == 0) {          if (e_offset0 == 0) {
349  #pragma omp parallel for private(i1,i2,k,node0)  #pragma omp parallel for
350              for (i2 = 0; i2 < local_NE2; i2++) {              for (index_t i2 = 0; i2 < local_NE2; i2++) {
351                  for (i1 = 0; i1 < local_NE1; i1++) {                  for (index_t i1 = 0; i1 < local_NE1; i1++) {
352                        const index_t k = 2 * (i1 + local_NE1 * i2) + faceNECount;
353                      index_t res, n0, n2, n4, n6;                      const index_t node0 = Nstride1 * (i1 + e_offset1)
354                      k = 2 * (i1 + local_NE1 * i2) + faceNECount;                                          + Nstride2 * (i2 + e_offset2);
355                      node0 = Nstride1 * N_PER_E * (i1 + e_offset1) + Nstride2 * N_PER_E * (i2 + e_offset2);                      const index_t res = 2 * (i1 + e_offset1)
356                      res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;                                        + NE1 * (i2 + e_offset2) + totalNECount;
357                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
358                      out->FaceElements->Tag[k] = LEFTTAG;                      out->FaceElements->Tag[k] = LEFTTAG;
359                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 370  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 361  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
361                      out->FaceElements->Tag[k + 1] = LEFTTAG;                      out->FaceElements->Tag[k + 1] = LEFTTAG;
362                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
363
364                      n0 = node0;                      const index_t n0 = node0;
365                      n2 = node0 + Nstride1;                      const index_t n2 = node0 + Nstride1;
366                      n4 = node0 + Nstride2;                      const index_t n4 = node0 + Nstride2;
367                      n6 = node0 + Nstride1 + Nstride2;                      const index_t n6 = node0 + Nstride1 + Nstride2;
368
369                      if ((global_adjustment + 0 + i1 + i2) % 2 == 0) {                      if ((global_adjustment + 0 + i1 + i2) % 2 == 0) {
370                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
# Line 384  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 375  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
375                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n6;
376                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n2;
377                      } else {                      } else {
378                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                          // this form is rotated around the 0,2,4,6 face
379                            // clockwise 90 degrees
380                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
381                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n4;
382                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n2;
# Line 399  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 390  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
390              faceNECount += 2 * local_NE1 * local_NE2;              faceNECount += 2 * local_NE1 * local_NE2;
391          }          }
392          totalNECount += 2 * NE1 * NE2;          totalNECount += 2 * NE1 * NE2;
393          /* **  elements on boundary 002 (x1=1): - Right */          // ** elements on boundary 002 (x1=1) - Right
394          if (local_NE0 + e_offset0 == NE0) {          if (local_NE0 + e_offset0 == NE0) {
395  #pragma omp parallel for private(i1,i2,k,node0)  #pragma omp parallel for
396              for (i2 = 0; i2 < local_NE2; i2++) {              for (index_t i2 = 0; i2 < local_NE2; i2++) {
397                  for (i1 = 0; i1 < local_NE1; i1++) {                  for (index_t i1 = 0; i1 < local_NE1; i1++) {
398                      index_t res, n1, n3, n5, n7;                      const index_t k = 2 * (i1 + local_NE1 * i2) + faceNECount;
399                      k = 2 * (i1 + local_NE1 * i2) + faceNECount;                      const index_t node0 = Nstride0 * (NE0 - 1)
400                                            + Nstride1 * (i1 + e_offset1)
401                      node0 =                                          + Nstride2 * (i2 + e_offset2);
402                          Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1) +                      const index_t res = 2 * (i1 + e_offset1)
403                          Nstride2 * N_PER_E * (i2 + e_offset2);                                        + NE1 * (i2 + e_offset2) + totalNECount;
res = 2 * ((i1 + e_offset1) + NE1 * (i2 + e_offset2)) + totalNECount;
404                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
405                      out->FaceElements->Tag[k] = RIGHTTAG;                      out->FaceElements->Tag[k] = RIGHTTAG;
406                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 418  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 408  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
408                      out->FaceElements->Tag[k + 1] = RIGHTTAG;                      out->FaceElements->Tag[k + 1] = RIGHTTAG;
409                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
410
411                      n1 = node0 + Nstride0;                      const index_t n1 = node0 + Nstride0;
412                      n3 = node0 + Nstride0 + Nstride1;                      const index_t n3 = node0 + Nstride0 + Nstride1;
413                      n5 = node0 + Nstride0 + Nstride2;                      const index_t n5 = node0 + Nstride0 + Nstride2;
414                      n7 = node0 + Nstride0 + Nstride1 + Nstride2;                      const index_t n7 = node0 + Nstride0 + Nstride1 + Nstride2;
415
416                      if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0) {                      if ((global_adjustment + local_NE0 - 1 + i1 + i2) % 2 == 0) {
417                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
418                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n3;
# Line 431  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 422  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
422                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
423                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n5;
424                      } else {                      } else {
425                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                          // this form is rotated around the 0,2,4,6 face
426                            // clockwise 90 degrees
427                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n1;
428                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7;                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n7;
429                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n5;
# Line 448  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 439  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
439          totalNECount += 2 * NE1 * NE2;          totalNECount += 2 * NE1 * NE2;
440      }      }
441      if (local_NE1 > 0) {      if (local_NE1 > 0) {
442          /* **  elements on boundary 010 (x2=0): -Front */          // ** elements on boundary 010 (x2=0) - Front
443          if (e_offset1 == 0) {          if (e_offset1 == 0) {
444  #pragma omp parallel for private(i0,i2,k,node0)  #pragma omp parallel for
445              for (i2 = 0; i2 < local_NE2; i2++) {              for (index_t i2 = 0; i2 < local_NE2; i2++) {
446                  for (i0 = 0; i0 < local_NE0; i0++) {                  for (index_t i0 = 0; i0 < local_NE0; i0++) {
447                      index_t res, n0, n1, n4, n5;                      const index_t k = 2 * (i0 + local_NE0 * i2) + faceNECount;
448                      k = 2 * (i0 + local_NE0 * i2) + faceNECount;                      const index_t node0 = Nstride0 * (i0 + e_offset0)
449                      node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride2 * N_PER_E * (i2 + e_offset2);                                          + Nstride2 * (i2 + e_offset2);
450                      res = 2 * ((i2 + e_offset2) + NE2 * (e_offset0 + i0)) + totalNECount;                      const index_t res = 2 * (i2 + e_offset2)
451                                          + NE2 * (e_offset0 + i0) + totalNECount;
452                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
453                      out->FaceElements->Tag[k] = FRONTTAG;                      out->FaceElements->Tag[k] = FRONTTAG;
454                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 464  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 456  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
456                      out->FaceElements->Tag[k + 1] = FRONTTAG;                      out->FaceElements->Tag[k + 1] = FRONTTAG;
457                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
458
459                      n0 = node0;                      const index_t n0 = node0;
460                      n1 = node0 + Nstride0;                      const index_t n1 = node0 + Nstride0;
461                      n4 = node0 + Nstride2;                      const index_t n4 = node0 + Nstride2;
462                      n5 = node0 + Nstride0 + Nstride2;                      const index_t n5 = node0 + Nstride0 + Nstride2;
463
464                      if ((global_adjustment + i0 + 0 + i2) % 2 == 0) {                      if ((global_adjustment + i0 + 0 + i2) % 2 == 0) {
465                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
# Line 479  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 471  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
471                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n4;
472
473                      } else {                      } else {
474                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                          // this form is rotated around the 0,2,4,6 face
475                            // clockwise 90 degrees
476                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n0;
477                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n1;
478                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4;                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n4;
# Line 495  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 487  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
487              faceNECount += 2 * local_NE0 * local_NE2;              faceNECount += 2 * local_NE0 * local_NE2;
488          }          }
489          totalNECount += 2 * NE0 * NE2;          totalNECount += 2 * NE0 * NE2;
490          /* **  elements on boundary 020 (x2=1): - Back */          // ** elements on boundary 020 (x2=1) - Back
491          if (local_NE1 + e_offset1 == NE1) {          if (local_NE1 + e_offset1 == NE1) {
492  #pragma omp parallel for private(i0,i2,k,node0)  #pragma omp parallel for
493              for (i2 = 0; i2 < local_NE2; i2++) {              for (index_t i2 = 0; i2 < local_NE2; i2++) {
494                  for (i0 = 0; i0 < local_NE0; i0++) {                  for (index_t i0 = 0; i0 < local_NE0; i0++) {
495                      index_t res, n2, n6, n7, n3;                      const index_t k = 2 * (i0 + local_NE0 * i2) + faceNECount;
496                      k = 2 * (i0 + local_NE0 * i2) + faceNECount;                      const index_t node0 = Nstride0 * (i0 + e_offset0)
497                      node0 =                                          + Nstride1 * (NE1 - 1)
498                          Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1) +                                          + Nstride2 * (i2 + e_offset2);
499                          Nstride2 * N_PER_E * (i2 + e_offset2);                      const index_t res = 2 * (i2 + e_offset2)
500                      res = 2 * ((i2 + e_offset2) + NE2 * (i0 + e_offset0)) + totalNECount;                                        + NE2 * (i0 + e_offset0) + totalNECount;
501                      out->FaceElements->Id[k] = res;                      out->FaceElements->Id[k] = res;
502                      out->FaceElements->Tag[k] = BACKTAG;                      out->FaceElements->Tag[k] = BACKTAG;
503                      out->FaceElements->Owner[k] = myRank;                      out->FaceElements->Owner[k] = myRank;
# Line 513  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 505  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
505                      out->FaceElements->Tag[k + 1] = BACKTAG;                      out->FaceElements->Tag[k + 1] = BACKTAG;
506                      out->FaceElements->Owner[k + 1] = myRank;                      out->FaceElements->Owner[k + 1] = myRank;
507
508                      n2 = node0 + Nstride1;                      const index_t n2 = node0 + Nstride1;
509                      n6 = node0 + Nstride1 + Nstride2;                      const index_t n6 = node0 + Nstride1 + Nstride2;
510                      n7 = node0 + Nstride0 + Nstride1 + Nstride2;                      const index_t n7 = node0 + Nstride0 + Nstride1 + Nstride2;
511                      n3 = node0 + Nstride0 + Nstride1;                      const index_t n3 = node0 + Nstride0 + Nstride1;
512
513                      if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0) {                      if ((global_adjustment + i0 + local_NE1 - 1 + i2) % 2 == 0) {
514                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
# Line 528  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 520  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
520                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;
521
522                      } else {                      } else {
523                          /* this form is rotated around the 0,2,4,6 face clockwise 90 degrees */                          // this form is rotated around the 0,2,4,6 face
524                            // clockwise 90 degrees
525                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(0, k, NN)] = n2;
526                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;                          out->FaceElements->Nodes[INDEX2(1, k, NN)] = n6;
527                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;                          out->FaceElements->Nodes[INDEX2(2, k, NN)] = n7;
# Line 536  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 529  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
529                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2;                          out->FaceElements->Nodes[INDEX2(0, k + 1, NN)] = n2;
530                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;                          out->FaceElements->Nodes[INDEX2(1, k + 1, NN)] = n7;
531                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;                          out->FaceElements->Nodes[INDEX2(2, k + 1, NN)] = n3;

532                      }                      }
533                  }                  }
534              }              }
# Line 544  Dudley_Mesh *Dudley_TriangularMesh_Tet4( Line 536  Dudley_Mesh *Dudley_TriangularMesh_Tet4(
536          }          }
537          totalNECount += 2 * NE0 * NE2;          totalNECount += 2 * NE0 * NE2;
538      }      }
// prepare mesh for further calculations
Dudley_Mesh_resolveNodeIds(out);
Dudley_Mesh_prepare(out, optimize);
539