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

Diff of /branches/trilinos_from_5897/dudley/src/Mesh_optimizeDOFLabeling.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: Mesh: optimizes the labeling of the DOFs on each processor */  
   
 /****************************************************************************/  
   
17  #include "Mesh.h"  #include "Mesh.h"
18  #include "IndexList.h"  #include "IndexList.h"
19    
# Line 27  Line 21 
21    
22  namespace dudley {  namespace dudley {
23    
24  void Dudley_Mesh_optimizeDOFLabeling(Dudley_Mesh* in, dim_t* distribution)  /// optimizes the labeling of the DOFs on each processor
25    void Mesh::optimizeDOFLabeling(const std::vector<index_t>& distribution)
26  {  {
27      if (in == NULL)      const int myRank = MPIInfo->rank;
28          return;      const int mpiSize = MPIInfo->size;
29      if (in->Nodes == NULL)      const index_t myFirstVertex = distribution[myRank];
30          return;      const index_t myLastVertex = distribution[myRank+1];
31        const dim_t myNumVertices = myLastVertex-myFirstVertex;
32      index_t myFirstVertex, myLastVertex, *newGlobalDOFID = NULL, firstVertex, lastVertex;      dim_t len = 0;
33      index_t k;      for (int p=0; p<mpiSize; ++p)
34      dim_t mpiSize, myNumVertices, len, p, i;          len=std::max(len, distribution[p+1]-distribution[p]);
     paso::Pattern_ptr pattern;  
     int myRank, current_rank;  
 #ifdef ESYS_MPI  
     int dest, source;  
     MPI_Status status;  
 #endif  
   
     myRank = in->MPIInfo->rank;  
     mpiSize = in->MPIInfo->size;  
     myFirstVertex = distribution[myRank];  
     myLastVertex = distribution[myRank + 1];  
     myNumVertices = myLastVertex - myFirstVertex;  
     len = 0;  
     for (p = 0; p < mpiSize; ++p)  
         len = std::max(len, distribution[p + 1] - distribution[p]);  
35    
36      boost::scoped_array<IndexList> index_list(new IndexList[myNumVertices]);      boost::scoped_array<IndexList> index_list(new IndexList[myNumVertices]);
37      newGlobalDOFID = new index_t[len];      std::vector<index_t> newGlobalDOFID(len);
38    
39      /* create the adjacency structure xadj and adjncy */      // create the adjacency structure xadj and adjncy
40  #pragma omp parallel private(i)  #pragma omp parallel
41      {      {
42          /*  insert contributions from element matrices into columns index index_list: */          // insert contributions from element matrices into columns index
43          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),          IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
44              myFirstVertex, myLastVertex, in->Elements,              myFirstVertex, myLastVertex, Elements,
45              in->Nodes->globalDegreesOfFreedom, in->Nodes->globalDegreesOfFreedom);              Nodes->globalDegreesOfFreedom);
46          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),          IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
47              myFirstVertex, myLastVertex, in->FaceElements,              myFirstVertex, myLastVertex, FaceElements,
48              in->Nodes->globalDegreesOfFreedom,              Nodes->globalDegreesOfFreedom);
49              in->Nodes->globalDegreesOfFreedom);          IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),
50          Dudley_IndexList_insertElementsWithRowRangeNoMainDiagonal(index_list.get(),              myFirstVertex, myLastVertex, Points,
51              myFirstVertex, myLastVertex, in->Points,              Nodes->globalDegreesOfFreedom);
             in->Nodes->globalDegreesOfFreedom,  
             in->Nodes->globalDegreesOfFreedom);  
52      }      }
53      /* create the local matrix pattern */      // create the local matrix pattern
54      pattern = paso::Pattern::fromIndexListArray(0, myNumVertices, index_list.get(),      paso::Pattern_ptr pattern = paso::Pattern::fromIndexListArray(0,
55          myFirstVertex, myLastVertex, -myFirstVertex);              myNumVertices, index_list.get(), myFirstVertex, myLastVertex,
56                -myFirstVertex);
57      pattern->reduceBandwidth(newGlobalDOFID);  
58        pattern->reduceBandwidth(&newGlobalDOFID[0]);
59      /* shift new labeling to create a global id */  
60  #pragma omp parallel for private(i)      // shift new labeling to create a global id
61      for (i = 0; i < myNumVertices; ++i)  #pragma omp parallel for
62        for (index_t i = 0; i < myNumVertices; ++i)
63          newGlobalDOFID[i] += myFirstVertex;          newGlobalDOFID[i] += myFirstVertex;
64    
65      /* distribute new labeling to other processors */      // distribute new labeling to other processors
66  #ifdef ESYS_MPI  #ifdef ESYS_MPI
67      dest = in->MPIInfo->mod_rank(myRank + 1);      const int dest = MPIInfo->mod_rank(myRank + 1);
68      source = in->MPIInfo->mod_rank(myRank - 1);      const int source = MPIInfo->mod_rank(myRank - 1);
69  #endif  #endif
70      current_rank = myRank;      int current_rank = myRank;
71      for (p = 0; p < mpiSize; ++p) {      for (int p = 0; p < mpiSize; ++p) {
72          firstVertex = distribution[current_rank];          const index_t firstVertex = distribution[current_rank];
73          lastVertex = distribution[current_rank + 1];          const index_t lastVertex = distribution[current_rank + 1];
74  #pragma omp parallel for private(i,k)  #pragma omp parallel for
75          for (i = 0; i < in->Nodes->numNodes; ++i) {          for (index_t i = 0; i < Nodes->getNumNodes(); ++i) {
76              k = in->Nodes->globalDegreesOfFreedom[i];              const index_t k = Nodes->globalDegreesOfFreedom[i];
77              if ((firstVertex <= k) && (k < lastVertex)) {              if (firstVertex <= k && k < lastVertex) {
78                  in->Nodes->globalDegreesOfFreedom[i] = newGlobalDOFID[k - firstVertex];                  Nodes->globalDegreesOfFreedom[i]=newGlobalDOFID[k-firstVertex];
79              }              }
80          }          }
81    
82          if (p < mpiSize - 1) { /* the final send can be skipped */          if (p < mpiSize - 1) { // the final send can be skipped
83  #ifdef ESYS_MPI  #ifdef ESYS_MPI
84              MPI_Sendrecv_replace(newGlobalDOFID, len, MPI_INT,              MPI_Status status;
85                                   dest, in->MPIInfo->counter(),              MPI_Sendrecv_replace(&newGlobalDOFID[0], len, MPI_DIM_T,
86                                   source, in->MPIInfo->counter(), in->MPIInfo->comm, &status);                                   dest, MPIInfo->counter(), source,
87              in->MPIInfo->incCounter();                                   MPIInfo->counter(), MPIInfo->comm, &status);
88                MPIInfo->incCounter();
89  #endif  #endif
90              current_rank = in->MPIInfo->mod_rank(current_rank - 1);              current_rank = MPIInfo->mod_rank(current_rank - 1);
91          }          }
92      }      }
     delete[] newGlobalDOFID;  
93  #if 0  #if 0
94      for (i = 0; i < in->Nodes->numNodes; ++i)      for (index_t i = 0; i < Nodes->getNumNodes(); ++i)
95          printf("%d ", in->Nodes->globalDegreesOfFreedom[i]);          printf("%d ", Nodes->globalDegreesOfFreedom[i]);
96      printf("\n");      printf("\n");
97  #endif  #endif
98  }  }

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

  ViewVC Help
Powered by ViewVC 1.1.26