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

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 triangular meshes by splitting rectangles */  /*   Dudley: generates triangular meshes by splitting rectangles */
20
21  /*   Generates a numElements[0] x numElements[1] x 2 mesh with first order elements (Tri3) in the rectangle */  /*   Generates a numElements[0] x numElements[1] x 2 mesh with first order elements (Tri3) in the rectangle */
22  /*   [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */  /*   [0,Length[0]] x [0,Length[1]]. order is the desired accuracy of the integration scheme. */
23
24  /************************************************************************************/  /****************************************************************************/

#define ESNEEDPYTHON
#include "esysUtils/first.h"
25
26  #include "TriangularMesh.h"  #include "TriangularMesh.h"
27
28  Dudley_Mesh *Dudley_TriangularMesh_Tri3(dim_t * numElements,  namespace dudley {
29                      double *Length, index_t order, index_t reduced_order, bool optimize, esysUtils::JMPI& mpi_info)
30    Dudley_Mesh *Dudley_TriangularMesh_Tri3(dim_t* numElements, double* Length,
31                     index_t order, index_t reduced_order, bool optimize,
32                     escript::JMPI& mpi_info)
33  {  {
34  #define N_PER_E 1  #define N_PER_E 1
35  #define DIM 2  #define DIM 2
# Line 39  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 39  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
39      index_t myRank;      index_t myRank;
40      Dudley_Mesh *out;      Dudley_Mesh *out;
41      char name[50];      char name[50];
42      const int LEFTTAG = 1;  /* boundary x1=0 */      const int LEFTTAG = 1;      /* boundary x1=0 */
43      const int RIGHTTAG = 2; /* boundary x1=1 */      const int RIGHTTAG = 2;     /* boundary x1=1 */
44      const int BOTTOMTAG = 10;   /* boundary x2=0 */      const int BOTTOMTAG = 10;   /* boundary x2=0 */
45      const int TOPTAG = 20;  /* boundary x2=1 */      const int TOPTAG = 20;      /* boundary x2=1 */
46
47  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
48      double time0 = Dudley_timer();      double time0 = Dudley_timer();
# Line 53  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 53  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
53
54      /* set up the global dimensions of the mesh */      /* set up the global dimensions of the mesh */
55
56      NE0 = MAX(1, numElements[0]);      NE0 = std::max(1, numElements[0]);
57      NE1 = MAX(1, numElements[1]);      NE1 = std::max(1, numElements[1]);
58      N0 = N_PER_E * NE0 + 1;      N0 = N_PER_E * NE0 + 1;
59      N1 = N_PER_E * NE1 + 1;      N1 = N_PER_E * NE1 + 1;
60
61      /* This code was originally copied from Finley's Rec4 constructor.      /* This code was originally copied from Finley's Rec4 constructor.
62         NE? refers to the number of rectangular elements in each direction.         NE? refers to the number of rectangular elements in each direction.
63         The number of nodes produced is the same but the number of non-face elements         The number of nodes produced is the same but the number of non-face
64         will double.         elements will double.
65       */       */
66
67      /*  allocate mesh: */      /*  allocate mesh: */
68      sprintf(name, "Triangular %d x %d (x 2) mesh", N0, N1);      sprintf(name, "Triangular %d x %d (x 2) mesh", N0, N1);
69      out = Dudley_Mesh_alloc(name, DIM, mpi_info);      out = Dudley_Mesh_alloc(name, DIM, mpi_info);
if (!Dudley_noError())
{
return NULL;
}
if (Dudley_noError())
{

Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));
Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Line2, mpi_info));
Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));
Nstride0 = 1;
Nstride1 = N0;
if (N1 == MAX(N0, N1))
{
local_NE0 = NE0;
e_offset0 = 0;
mpi_info->split(NE1, &local_NE1, &e_offset1);
}
else
{
mpi_info->split(NE0, &local_NE0, &e_offset0);
local_NE1 = NE1;
e_offset1 = 0;
}
offset0 = e_offset0 * N_PER_E;
offset1 = e_offset1 * N_PER_E;
local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;
local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;

/* get the number of surface elements */

NFaceElements = 0;
if (local_NE0 > 0)
{
NDOF0 = N0;
if (e_offset0 == 0)
NFaceElements += local_NE1;
if (local_NE0 + e_offset0 == NE0)
NFaceElements += local_NE1;
}
else
{
NDOF0 = N0 - 1;
}
if (local_NE1 > 0)
{
NDOF1 = N1;
if (e_offset1 == 0)
NFaceElements += local_NE0;
if (local_NE1 + e_offset1 == NE1)
NFaceElements += local_NE0;
}
else
{
NDOF1 = N1 - 1;
}

/*  allocate tables: */

Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1);

/* This code was originally copied from Finley's rec4 generator
We double these numbers because each "rectangle" will be split into
two triangles. So the number of nodes is the same but the
number of elements will double */
Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * 2);
Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);
70
71      }      Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));
72      if (Dudley_noError())      Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Line2, mpi_info));
73      {      Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));
74      dim_t NN;      Nstride0 = 1;
75      index_t global_adjustment;      Nstride1 = N0;
76      /* create nodes */      if (N1 == std::max(N0, N1)) {
77            local_NE0 = NE0;
78            e_offset0 = 0;
79            mpi_info->split(NE1, &local_NE1, &e_offset1);
80        } else {
81            mpi_info->split(NE0, &local_NE0, &e_offset0);
82            local_NE1 = NE1;
83            e_offset1 = 0;
84        }
85        offset0 = e_offset0 * N_PER_E;
86        offset1 = e_offset1 * N_PER_E;
87        local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;
88        local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;
89
90        /* get the number of surface elements */
91
92        NFaceElements = 0;
93        if (local_NE0 > 0) {
94            NDOF0 = N0;
95            if (e_offset0 == 0)
96                NFaceElements += local_NE1;
97            if (local_NE0 + e_offset0 == NE0)
98                NFaceElements += local_NE1;
99        } else {
100            NDOF0 = N0 - 1;
101        }
102        if (local_NE1 > 0) {
103            NDOF1 = N1;
104            if (e_offset1 == 0)
105                NFaceElements += local_NE0;
106            if (local_NE1 + e_offset1 == NE1)
107                NFaceElements += local_NE0;
108        } else {
109            NDOF1 = N1 - 1;
110        }
111
112        Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1);
113
114        /* This code was originally copied from Finley's rec4 generator
115           We double these numbers because each "rectangle" will be split into
116           two triangles. So the number of nodes is the same but the
117           number of elements will double */
118        Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * 2);
119        Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);
120
121        dim_t NN;
123        /* create nodes */
124  #pragma omp parallel for private(i0,i1)  #pragma omp parallel for private(i0,i1)
125      for (i1 = 0; i1 < local_N1; i1++)      for (i1 = 0; i1 < local_N1; i1++) {
126      {          for (i0 = 0; i0 < local_N0; i0++) {
127          for (i0 = 0; i0 < local_N0; i0++)              dim_t k = i0 + local_N0 * i1;
128          {              dim_t global_i0 = i0 + offset0;
129          dim_t k = i0 + local_N0 * i1;              dim_t global_i1 = i1 + offset1;
130          dim_t global_i0 = i0 + offset0;              out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0];
131          dim_t global_i1 = i1 + offset1;              out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];
132          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;
133          out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];              out->Nodes->Tag[k] = 0;
134          out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;              out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0) + Nstride1 * (global_i1 % NDOF1);
135          out->Nodes->Tag[k] = 0;          }
136          out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0) + Nstride1 * (global_i1 % NDOF1);      }
137          }      /*   set the elements: */
138      }      NN = out->Elements->numNodes;
139      /*   set the elements: */      global_adjustment = (offset0 + offset1) % 2;
NN = out->Elements->numNodes;
global_adjustment = (offset0 + offset1) % 2;
140  #pragma omp parallel for private(i0,i1)  #pragma omp parallel for private(i0,i1)
141      for (i1 = 0; i1 < local_NE1; i1++)      for (i1 = 0; i1 < local_NE1; i1++) {
142      {          for (i0 = 0; i0 < local_NE0; i0++) {
143          for (i0 = 0; i0 < local_NE0; i0++)              index_t a, b, c, d;
144          {              /* we will split this "rectangle" into two triangles */
145          index_t a, b, c, d;              dim_t k = 2 * (i0 + local_NE0 * i1);
146          /* we will split this "rectangle" into two triangles */              index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);
147          dim_t k = 2 * (i0 + local_NE0 * i1);
148          index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);              out->Elements->Id[k] = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1));
149                out->Elements->Tag[k] = 0;
150          out->Elements->Id[k] = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1));              out->Elements->Owner[k] = myRank;
151          out->Elements->Tag[k] = 0;              out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;
152          out->Elements->Owner[k] = myRank;              out->Elements->Tag[k + 1] = 0;
153          out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;              out->Elements->Owner[k + 1] = myRank;
154          out->Elements->Tag[k + 1] = 0;
155          out->Elements->Owner[k + 1] = myRank;              /* a,b,c,d gives the nodes in the rectangle in clockwise order */
156                a = node0; b = node0 + Nstride0; c = node0 + Nstride1 + Nstride0; d = node0 + Nstride1;
157          /* a,b,c,d gives the nodes in the rectangle in clockwise order */              /* For a little bit of variety  */
158          a = node0; b = node0 + Nstride0; c = node0 + Nstride1 + Nstride0; d = node0 + Nstride1;              if ((global_adjustment + node0) % 2) {
159          /* For a little bit of variety  */                  out->Elements->Nodes[INDEX2(0, k, NN)] = a;
160          if ((global_adjustment + node0) % 2)                  out->Elements->Nodes[INDEX2(1, k, NN)] = b;
161          {                  out->Elements->Nodes[INDEX2(2, k, NN)] = d;
162              out->Elements->Nodes[INDEX2(0, k, NN)] = a;                  out->Elements->Nodes[INDEX2(0, k + 1, NN)] = b;
163              out->Elements->Nodes[INDEX2(1, k, NN)] = b;                  out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
164              out->Elements->Nodes[INDEX2(2, k, NN)] = d;                  out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
165              out->Elements->Nodes[INDEX2(0, k + 1, NN)] = b;              } else {
166              out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;                  out->Elements->Nodes[INDEX2(0, k, NN)] = a;
167              out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;                  out->Elements->Nodes[INDEX2(1, k, NN)] = b;
168          }                  out->Elements->Nodes[INDEX2(2, k, NN)] = c;
169          else                  out->Elements->Nodes[INDEX2(0, k + 1, NN)] = a;
170          {                  out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;
171              out->Elements->Nodes[INDEX2(0, k, NN)] = a;                  out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;
172              out->Elements->Nodes[INDEX2(1, k, NN)] = b;              }
173              out->Elements->Nodes[INDEX2(2, k, NN)] = c;          }
174              out->Elements->Nodes[INDEX2(0, k + 1, NN)] = a;      }
175              out->Elements->Nodes[INDEX2(1, k + 1, NN)] = c;      /* face elements */
176              out->Elements->Nodes[INDEX2(2, k + 1, NN)] = d;      NN = out->FaceElements->numNodes;
177          }      totalNECount = 2 * NE0 * NE1;   /* because we have split the rectangles */
178          }      faceNECount = 0;
179      }      if (local_NE0 > 0) {
180      /* face elements */          /* **  elements on boundary 001 (x1=0): */
NN = out->FaceElements->numNodes;
totalNECount = 2 * NE0 * NE1;   /* because we have split the rectangles */
faceNECount = 0;
if (local_NE0 > 0)
{
/* **  elements on boundary 001 (x1=0): */
181
182          if (e_offset0 == 0)          if (e_offset0 == 0) {
{
183  #pragma omp parallel for private(i1)  #pragma omp parallel for private(i1)
184          for (i1 = 0; i1 < local_NE1; i1++)              for (i1 = 0; i1 < local_NE1; i1++) {
185          {                  dim_t k = i1 + faceNECount;
186                    index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);
187              dim_t k = i1 + faceNECount;
188              index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);                  out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;
189                    out->FaceElements->Tag[k] = LEFTTAG;
190              out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;                  out->FaceElements->Owner[k] = myRank;
191              out->FaceElements->Tag[k] = LEFTTAG;                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1;
192              out->FaceElements->Owner[k] = myRank;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0;
193              out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1;              }
194              out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0;              faceNECount += local_NE1;
195          }          }
196          faceNECount += local_NE1;          totalNECount += NE1;
197          }          /* **  elements on boundary 002 (x1=1): */
198          totalNECount += NE1;          if (local_NE0 + e_offset0 == NE0) {
/* **  elements on boundary 002 (x1=1): */
if (local_NE0 + e_offset0 == NE0)
{
199  #pragma omp parallel for private(i1)  #pragma omp parallel for private(i1)
200          for (i1 = 0; i1 < local_NE1; i1++)              for (i1 = 0; i1 < local_NE1; i1++) {
201          {                  dim_t k = i1 + faceNECount;
202              dim_t k = i1 + faceNECount;                  index_t node0 = Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1);
203              index_t node0 = Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1);
204                    out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;
205              out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;                  out->FaceElements->Tag[k] = RIGHTTAG;
206              out->FaceElements->Tag[k] = RIGHTTAG;                  out->FaceElements->Owner[k] = myRank;
207              out->FaceElements->Owner[k] = myRank;                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride0;
208              out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride0;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;
209              out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1 + Nstride0;              }
210          }              faceNECount += local_NE1;
211          faceNECount += local_NE1;          }
212          }          totalNECount += NE1;
213          totalNECount += NE1;      }
214      }      if (local_NE1 > 0) {
215      if (local_NE1 > 0)          /* **  elements on boundary 010 (x2=0): */
216      {          if (e_offset1 == 0) {
/* **  elements on boundary 010 (x2=0): */
if (e_offset1 == 0)
{
217  #pragma omp parallel for private(i0)  #pragma omp parallel for private(i0)
218          for (i0 = 0; i0 < local_NE0; i0++)              for (i0 = 0; i0 < local_NE0; i0++) {
219          {                  dim_t k = i0 + faceNECount;
220              dim_t k = i0 + faceNECount;                  index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);
221              index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);
222                    out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;
223              out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;                  out->FaceElements->Tag[k] = BOTTOMTAG;
224              out->FaceElements->Tag[k] = BOTTOMTAG;                  out->FaceElements->Owner[k] = myRank;
225              out->FaceElements->Owner[k] = myRank;
226                    out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;
227              out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;
228              out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;              }
229          }              faceNECount += local_NE0;
230          faceNECount += local_NE0;          }
231          }          totalNECount += NE0;
232          totalNECount += NE0;          /* **  elements on boundary 020 (x2=1): */
233          /* **  elements on boundary 020 (x2=1): */          if (local_NE1 + e_offset1 == NE1) {
if (local_NE1 + e_offset1 == NE1)
{
234  #pragma omp parallel for private(i0)  #pragma omp parallel for private(i0)
235          for (i0 = 0; i0 < local_NE0; i0++)              for (i0 = 0; i0 < local_NE0; i0++) {
236          {                  dim_t k = i0 + faceNECount;
237              dim_t k = i0 + faceNECount;                  index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1);
238              index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1);
239                    out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;
240              out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;                  out->FaceElements->Tag[k] = TOPTAG;
241              out->FaceElements->Tag[k] = TOPTAG;                  out->FaceElements->Owner[k] = myRank;
out->FaceElements->Owner[k] = myRank;
242
243              out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
244              out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;
245  /*printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],  /*printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],
246  INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]); */  INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]); */
247          }              }
248          faceNECount += local_NE0;              faceNECount += local_NE0;
249          }          }
250          totalNECount += NE0;          totalNECount += NE0;
251      }      }
252      }      /* add tag names */
253      if (Dudley_noError())      Dudley_Mesh_addTagMap(out, "top", TOPTAG);
}
257      /* prepare mesh for further calculations: */      /* prepare mesh for further calculations: */
258      if (Dudley_noError())      Dudley_Mesh_resolveNodeIds(out);
259      {      Dudley_Mesh_prepare(out, optimize);
Dudley_Mesh_resolveNodeIds(out);
}
if (Dudley_noError())
{
Dudley_Mesh_prepare(out, optimize);
}
260      return out;      return out;
261  }  }
262
263    } // namespace dudley
264

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