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

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

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

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 triangular meshes by splitting rectangles */  
   
 /*   Generates a numElements[0] x numElements[1] x 2 mesh with first order elements (Tri3) in the rectangle */  
 /*   [0,Length[0]] x [0,Length[1]]. 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  namespace dudley {  namespace dudley {
22    
23  Dudley_Mesh *Dudley_TriangularMesh_Tri3(dim_t* numElements, double* Length,  Mesh* TriangularMesh_Tri3(const dim_t* numElements, const double* length,
24                   index_t order, index_t reduced_order, bool optimize,                            bool optimize, escript::JMPI mpiInfo)
                  escript::JMPI& mpi_info)  
25  {  {
26  #define N_PER_E 1      const int DIM = 2;
27  #define DIM 2      const int LEFTTAG = 1;    // boundary x1=0
28      dim_t N0, N1, NE0, NE1, i0, i1, Nstride0 = 0, Nstride1 = 0, local_NE0, local_NE1, local_N0 = 0, local_N1 = 0;      const int RIGHTTAG = 2;   // boundary x1=1
29      index_t offset0 = 0, offset1 = 0, e_offset0 = 0, e_offset1 = 0;      const int BOTTOMTAG = 10; // boundary x2=0
30      dim_t totalNECount, faceNECount, NDOF0 = 0, NDOF1 = 0, NFaceElements;      const int TOPTAG = 20;    // boundary x2=1
     index_t myRank;  
     Dudley_Mesh *out;  
     char name[50];  
     const int LEFTTAG = 1;      /* boundary x1=0 */  
     const int RIGHTTAG = 2;     /* boundary x1=1 */  
     const int BOTTOMTAG = 10;   /* boundary x2=0 */  
     const int TOPTAG = 20;      /* boundary x2=1 */  
31    
32  #ifdef Dudley_TRACE  #ifdef Dudley_TRACE
33      double time0 = Dudley_timer();      double time0 = Dudley_timer();
34  #endif  #endif
35    
36      /* get MPI information */      const int myRank = mpiInfo->rank;
     myRank = mpi_info->rank;  
37    
38      /* set up the global dimensions of the mesh */      // set up the global dimensions of the mesh
39        const dim_t NE0 = std::max(1, numElements[0]);
40      NE0 = std::max(1, numElements[0]);      const dim_t NE1 = std::max(1, numElements[1]);
41      NE1 = std::max(1, numElements[1]);      const dim_t N0 = NE0 + 1;
42      N0 = N_PER_E * NE0 + 1;      const dim_t N1 = NE1 + 1;
43      N1 = N_PER_E * NE1 + 1;  
44        // This code was originally copied from Finley's Rec4 constructor.
45      /* This code was originally copied from Finley's Rec4 constructor.      // NE? refers to the number of rectangular elements in each direction.
46         NE? refers to the number of rectangular elements in each direction.      // The number of nodes produced is the same but the number of non-face
47         The number of nodes produced is the same but the number of non-face      // elements will double since each "rectangle" is split into two triangles.
48         elements will double.  
49       */      // allocate mesh
50        std::stringstream name;
51      /*  allocate mesh: */      name << "Triangular " << N0 << " x " << N1 << " (x 2) mesh";
52      sprintf(name, "Triangular %d x %d (x 2) mesh", N0, N1);      Mesh* out = new Mesh(name.str(), DIM, mpiInfo);
53      out = Dudley_Mesh_alloc(name, DIM, mpi_info);  
54        out->setPoints(new ElementFile(Dudley_Point1, mpiInfo));
55      Dudley_Mesh_setPoints(out, Dudley_ElementFile_alloc(Dudley_Point1, mpi_info));      out->setFaceElements(new ElementFile(Dudley_Line2, mpiInfo));
56      Dudley_Mesh_setFaceElements(out, Dudley_ElementFile_alloc(Dudley_Line2, mpi_info));      out->setElements(new ElementFile(Dudley_Tri3, mpiInfo));
57      Dudley_Mesh_setElements(out, Dudley_ElementFile_alloc(Dudley_Tri3, mpi_info));  
58      Nstride0 = 1;      const dim_t Nstride0 = 1;
59      Nstride1 = N0;      const dim_t Nstride1 = N0;
60        dim_t local_NE0, local_NE1;
61        index_t e_offset0 = 0, e_offset1 = 0;
62      if (N1 == std::max(N0, N1)) {      if (N1 == std::max(N0, N1)) {
63          local_NE0 = NE0;          local_NE0 = NE0;
64          e_offset0 = 0;          e_offset0 = 0;
65          mpi_info->split(NE1, &local_NE1, &e_offset1);          mpiInfo->split(NE1, &local_NE1, &e_offset1);
66      } else {      } else {
67          mpi_info->split(NE0, &local_NE0, &e_offset0);          mpiInfo->split(NE0, &local_NE0, &e_offset0);
68          local_NE1 = NE1;          local_NE1 = NE1;
69          e_offset1 = 0;          e_offset1 = 0;
70      }      }
71      offset0 = e_offset0 * N_PER_E;      const index_t offset0 = e_offset0;
72      offset1 = e_offset1 * N_PER_E;      const index_t offset1 = e_offset1;
73      local_N0 = local_NE0 > 0 ? local_NE0 * N_PER_E + 1 : 0;      const dim_t local_N0 = local_NE0 > 0 ? local_NE0 + 1 : 0;
74      local_N1 = local_NE1 > 0 ? local_NE1 * N_PER_E + 1 : 0;      const dim_t local_N1 = local_NE1 > 0 ? local_NE1 + 1 : 0;
75    
76      /* get the number of surface elements */      // get the number of surface elements
77        dim_t NFaceElements = 0;
78      NFaceElements = 0;      dim_t NDOF0, NDOF1;
79      if (local_NE0 > 0) {      if (local_NE0 > 0) {
80          NDOF0 = N0;          NDOF0 = N0;
81          if (e_offset0 == 0)          if (e_offset0 == 0)
# Line 109  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 95  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
95          NDOF1 = N1 - 1;          NDOF1 = N1 - 1;
96      }      }
97    
98      Dudley_NodeFile_allocTable(out->Nodes, local_N0 * local_N1);      out->Nodes->allocTable(local_N0 * local_N1);
99        out->Elements->allocTable(local_NE0 * local_NE1 * 2);
100      /* This code was originally copied from Finley's rec4 generator      out->FaceElements->allocTable(NFaceElements);
101         We double these numbers because each "rectangle" will be split into  
102         two triangles. So the number of nodes is the same but the      // create nodes
103         number of elements will double */  #pragma omp parallel for
104      Dudley_ElementFile_allocTable(out->Elements, local_NE0 * local_NE1 * 2);      for (index_t i1 = 0; i1 < local_N1; i1++) {
105      Dudley_ElementFile_allocTable(out->FaceElements, NFaceElements);          for (index_t i0 = 0; i0 < local_N0; i0++) {
106                const dim_t k = i0 + local_N0 * i1;
107      dim_t NN;              const dim_t global_i0 = i0 + offset0;
108      index_t global_adjustment;              const dim_t global_i1 = i1 + offset1;
109      /* create nodes */              out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (real_t)global_i0 / (real_t)(N0 - 1) * length[0];
110  #pragma omp parallel for private(i0,i1)              out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (real_t)global_i1 / (real_t)(N1 - 1) * length[1];
     for (i1 = 0; i1 < local_N1; i1++) {  
         for (i0 = 0; i0 < local_N0; i0++) {  
             dim_t k = i0 + local_N0 * i1;  
             dim_t global_i0 = i0 + offset0;  
             dim_t global_i1 = i1 + offset1;  
             out->Nodes->Coordinates[INDEX2(0, k, DIM)] = (double)global_i0 / (double)(N0 - 1) * Length[0];  
             out->Nodes->Coordinates[INDEX2(1, k, DIM)] = (double)global_i1 / (double)(N1 - 1) * Length[1];  
111              out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;              out->Nodes->Id[k] = Nstride0 * global_i0 + Nstride1 * global_i1;
112              out->Nodes->Tag[k] = 0;              out->Nodes->Tag[k] = 0;
113              out->Nodes->globalDegreesOfFreedom[k] = Nstride0 * (global_i0 % NDOF0) + Nstride1 * (global_i1 % NDOF1);              out->Nodes->globalDegreesOfFreedom[k] =
114                                                    Nstride0 * (global_i0 % NDOF0)
115                                                  + Nstride1 * (global_i1 % NDOF1);
116          }          }
117      }      }
     /*   set the elements: */  
     NN = out->Elements->numNodes;  
     global_adjustment = (offset0 + offset1) % 2;  
 #pragma omp parallel for private(i0,i1)  
     for (i1 = 0; i1 < local_NE1; i1++) {  
         for (i0 = 0; i0 < local_NE0; i0++) {  
             index_t a, b, c, d;  
             /* we will split this "rectangle" into two triangles */  
             dim_t k = 2 * (i0 + local_NE0 * i1);  
             index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (i1 + e_offset1);  
118    
119              out->Elements->Id[k] = 2 * ((i0 + e_offset0) + NE0 * (i1 + e_offset1));      // set the elements
120        dim_t NN = out->Elements->numNodes;
121        const index_t global_adjustment = (offset0 + offset1) % 2;
122    
123    #pragma omp parallel for
124        for (index_t i1 = 0; i1 < local_NE1; i1++) {
125            for (index_t i0 = 0; i0 < local_NE0; i0++) {
126                // we will split this "rectangle" into two triangles
127                const dim_t k = 2 * (i0 + local_NE0 * i1);
128                const index_t node0 = Nstride0 * (i0 + e_offset0)
129                                    + Nstride1 * (i1 + e_offset1);
130    
131                out->Elements->Id[k] = 2 * (i0 + e_offset0) + NE0*(i1 + e_offset1);
132              out->Elements->Tag[k] = 0;              out->Elements->Tag[k] = 0;
133              out->Elements->Owner[k] = myRank;              out->Elements->Owner[k] = myRank;
134              out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;              out->Elements->Id[k + 1] = out->Elements->Id[k] + 1;
135              out->Elements->Tag[k + 1] = 0;              out->Elements->Tag[k + 1] = 0;
136              out->Elements->Owner[k + 1] = myRank;              out->Elements->Owner[k + 1] = myRank;
137    
138              /* a,b,c,d gives the nodes in the rectangle in clockwise order */              // a,b,c,d gives the nodes of the rectangle in clockwise order
139              a = node0; b = node0 + Nstride0; c = node0 + Nstride1 + Nstride0; d = node0 + Nstride1;              const index_t a = node0;
140              /* For a little bit of variety  */              const index_t b = node0 + Nstride0;
141                const index_t c = node0 + Nstride1 + Nstride0;
142                const index_t d = node0 + Nstride1;
143                // For a little bit of variety
144              if ((global_adjustment + node0) % 2) {              if ((global_adjustment + node0) % 2) {
145                  out->Elements->Nodes[INDEX2(0, k, NN)] = a;                  out->Elements->Nodes[INDEX2(0, k, NN)] = a;
146                  out->Elements->Nodes[INDEX2(1, k, NN)] = b;                  out->Elements->Nodes[INDEX2(1, k, NN)] = b;
# Line 172  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 158  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
158              }              }
159          }          }
160      }      }
161      /* face elements */  
162        // face elements
163      NN = out->FaceElements->numNodes;      NN = out->FaceElements->numNodes;
164      totalNECount = 2 * NE0 * NE1;   /* because we have split the rectangles */      dim_t totalNECount = 2 * NE0 * NE1; // because we have split the rectangles
165      faceNECount = 0;      dim_t faceNECount = 0;
166      if (local_NE0 > 0) {      if (local_NE0 > 0) {
167          /* **  elements on boundary 001 (x1=0): */          // ** elements on boundary 001 (x1=0)
   
168          if (e_offset0 == 0) {          if (e_offset0 == 0) {
169  #pragma omp parallel for private(i1)  #pragma omp parallel for
170              for (i1 = 0; i1 < local_NE1; i1++) {              for (index_t i1 = 0; i1 < local_NE1; i1++) {
171                  dim_t k = i1 + faceNECount;                  const dim_t k = i1 + faceNECount;
172                  index_t node0 = Nstride1 * N_PER_E * (i1 + e_offset1);                  const index_t node0 = Nstride1 * (i1 + e_offset1);
   
173                  out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;                  out->FaceElements->Id[k] = i1 + e_offset1 + totalNECount;
174                  out->FaceElements->Tag[k] = LEFTTAG;                  out->FaceElements->Tag[k] = LEFTTAG;
175                  out->FaceElements->Owner[k] = myRank;                  out->FaceElements->Owner[k] = myRank;
# Line 194  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 179  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
179              faceNECount += local_NE1;              faceNECount += local_NE1;
180          }          }
181          totalNECount += NE1;          totalNECount += NE1;
182          /* **  elements on boundary 002 (x1=1): */          // ** elements on boundary 002 (x1=1)
183          if (local_NE0 + e_offset0 == NE0) {          if (local_NE0 + e_offset0 == NE0) {
184  #pragma omp parallel for private(i1)  #pragma omp parallel for
185              for (i1 = 0; i1 < local_NE1; i1++) {              for (index_t i1 = 0; i1 < local_NE1; i1++) {
186                  dim_t k = i1 + faceNECount;                  const dim_t k = i1 + faceNECount;
187                  index_t node0 = Nstride0 * N_PER_E * (NE0 - 1) + Nstride1 * N_PER_E * (i1 + e_offset1);                  const index_t node0 = Nstride0 * (NE0 - 1)
188                                        + Nstride1 * (i1 + e_offset1);
189    
190                  out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;                  out->FaceElements->Id[k] = (i1 + e_offset1) + totalNECount;
191                  out->FaceElements->Tag[k] = RIGHTTAG;                  out->FaceElements->Tag[k] = RIGHTTAG;
# Line 212  Dudley_Mesh *Dudley_TriangularMesh_Tri3( Line 198  Dudley_Mesh *Dudley_TriangularMesh_Tri3(
198          totalNECount += NE1;          totalNECount += NE1;
199      }      }
200      if (local_NE1 > 0) {      if (local_NE1 > 0) {
201          /* **  elements on boundary 010 (x2=0): */          // ** elements on boundary 010 (x2=0)
202          if (e_offset1 == 0) {          if (e_offset1 == 0) {
203  #pragma omp parallel for private(i0)  #pragma omp parallel for
204              for (i0 = 0; i0 < local_NE0; i0++) {              for (index_t i0 = 0; i0 < local_NE0; i0++) {
205                  dim_t k = i0 + faceNECount;                  const dim_t k = i0 + faceNECount;
206                  index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0);                  const index_t node0 = Nstride0 * (i0 + e_offset0);
   
207                  out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;                  out->FaceElements->Id[k] = e_offset0 + i0 + totalNECount;
208                  out->FaceElements->Tag[k] = BOTTOMTAG;                  out->FaceElements->Tag[k] = BOTTOMTAG;
209                  out->FaceElements->Owner[k] = myRank;                  out->FaceElements->Owner[k] = myRank;
   
210                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0;
211                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride0;
212              }              }
213              faceNECount += local_NE0;              faceNECount += local_NE0;
214          }          }
215          totalNECount += NE0;          totalNECount += NE0;
216          /* **  elements on boundary 020 (x2=1): */          // ** elements on boundary 020 (x2=1)
217          if (local_NE1 + e_offset1 == NE1) {          if (local_NE1 + e_offset1 == NE1) {
218  #pragma omp parallel for private(i0)  #pragma omp parallel for
219              for (i0 = 0; i0 < local_NE0; i0++) {              for (index_t i0 = 0; i0 < local_NE0; i0++) {
220                  dim_t k = i0 + faceNECount;                  const dim_t k = i0 + faceNECount;
221                  index_t node0 = Nstride0 * N_PER_E * (i0 + e_offset0) + Nstride1 * N_PER_E * (NE1 - 1);                  const index_t node0 = Nstride0 * (i0 + e_offset0)
222                                        + Nstride1 * (NE1 - 1);
223    
224                  out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;                  out->FaceElements->Id[k] = i0 + e_offset0 + totalNECount;
225                  out->FaceElements->Tag[k] = TOPTAG;                  out->FaceElements->Tag[k] = TOPTAG;
226                  out->FaceElements->Owner[k] = myRank;                  out->FaceElements->Owner[k] = myRank;
   
227                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;                  out->FaceElements->Nodes[INDEX2(0, k, NN)] = node0 + Nstride1 + Nstride0;
228                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;                  out->FaceElements->Nodes[INDEX2(1, k, NN)] = node0 + Nstride1;
 /*printf("E=%d: %d=%d %d=%d\n",k,INDEX2(0,k,NN),out->FaceElements->Nodes[INDEX2(0,k,NN)],  
 INDEX2(1,k,NN),out->FaceElements->Nodes[INDEX2(1,k,NN)]); */  
229              }              }
230              faceNECount += local_NE0;              faceNECount += local_NE0;
231          }          }
232          totalNECount += NE0;          totalNECount += NE0;
233      }      }
234      /* add tag names */  
235      Dudley_Mesh_addTagMap(out, "top", TOPTAG);      // add tag names
236      Dudley_Mesh_addTagMap(out, "bottom", BOTTOMTAG);      out->addTagMap("top", TOPTAG);
237      Dudley_Mesh_addTagMap(out, "left", LEFTTAG);      out->addTagMap("bottom", BOTTOMTAG);
238      Dudley_Mesh_addTagMap(out, "right", RIGHTTAG);      out->addTagMap("left", LEFTTAG);
239      /* prepare mesh for further calculations: */      out->addTagMap("right", RIGHTTAG);
240      Dudley_Mesh_resolveNodeIds(out);  
241      Dudley_Mesh_prepare(out, optimize);      // prepare mesh for further calculations
242        out->resolveNodeIds();
243        out->prepare(optimize);
244      return out;      return out;
245  }  }
246    

Legend:
Removed from v.6078  
changed lines
  Added in v.6079

  ViewVC Help
Powered by ViewVC 1.1.26