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

Diff of /branches/trilinos_from_5897/dudley/src/ElementFile.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    
17  /****************************************************************************/  #include "ElementFile.h"
18    #include "ShapeTable.h"
19    
20  /*   Dudley: ElementFile */  namespace dudley {
21    
22  /*   allocates an element file to hold elements of type id and with integration order order. */  ElementFile::ElementFile(ElementTypeId type, escript::JMPI mpiInfo) :
23  /*   use Dudley_Mesh_allocElementTable to allocate the element table (Id,Nodes,Tag,Owner). */      MPIInfo(mpiInfo),
24        numElements(0),
25        Id(NULL),
26        Tag(NULL),
27        Owner(NULL),
28        Nodes(NULL),
29        Color(NULL),
30        minColor(0),
31        maxColor(-1),
32        etype(type)
33    {
34        jacobians = new ElementFile_Jacobians();
35        jacobians_reducedQ = new ElementFile_Jacobians();
36    
37        numDim = Dims[type];
38        numNodes = numDim + 1;
39        numLocalDim = localDims[type];
40        numShapes = numLocalDim + 1;
41        ename = getElementName(type);
42    }
43    
44    ElementFile::~ElementFile()
45    {
46        freeTable();
47        delete jacobians;
48        delete jacobians_reducedQ;
49    }
50    
51    void ElementFile::allocTable(dim_t NE)
52    {
53        if (numElements > 0)
54            freeTable();
55    
56        numElements = NE;
57        Owner = new int[numElements];
58        Id = new index_t[numElements];
59        Nodes = new index_t[numElements * numNodes];
60        Tag = new int[numElements];
61        Color = new index_t[numElements];
62    
63        // this initialization makes sure that data are located on the right
64        // processor
65    #pragma omp parallel for
66        for (index_t e = 0; e < numElements; e++) {
67            for (int i = 0; i < numNodes; i++)
68                Nodes[INDEX2(i, e, numNodes)] = -1;
69            Owner[e] = -1;
70            Id[e] = -1;
71            Tag[e] = -1;
72            Color[e] = -1;
73        }
74        maxColor = -1;
75        minColor = 0;
76    }
77    
78  /****************************************************************************/  void ElementFile::freeTable()
79    {
80        delete[] Owner;
81        delete[] Id;
82        delete[] Nodes;
83        delete[] Tag;
84        delete[] Color;
85        tagsInUse.clear();
86        numElements = 0;
87        maxColor = -1;
88        minColor = 0;
89    }
90    
91  #include "ElementFile.h"  void ElementFile::copyTable(index_t offset, index_t nodeOffset,
92  #include "ShapeTable.h"                              index_t idOffset, const ElementFile* in)
93    {
94        const int NN_in = in->numNodes;
95        if (NN_in > numNodes) {
96            throw DudleyException("ElementFile::copyTable: dimensions of element files don't match.");
97        }
98    
99  namespace dudley {      if (MPIInfo->comm != in->MPIInfo->comm) {
100            throw DudleyException("ElementFile::copyTable: MPI communicators of element files don't match.");
101        }
102    
103  Dudley_ElementFile *Dudley_ElementFile_alloc(Dudley_ElementTypeId etype, escript::JMPI& MPIInfo)  #pragma omp parallel for
104        for (index_t n = 0; n < in->numElements; n++) {
105            Owner[offset + n] = in->Owner[n];
106            Id[offset + n] = in->Id[n] + idOffset;
107            Tag[offset + n] = in->Tag[n];
108            for (int i = 0; i < numNodes; i++)
109                Nodes[INDEX2(i, offset + n, numNodes)] =
110                                in->Nodes[INDEX2(i, n, NN_in)] + nodeOffset;
111        }
112    }
113    
114    void ElementFile::gather(const index_t* index, const ElementFile* in)
115  {  {
116      Dudley_ElementFile* out = new Dudley_ElementFile;      const int NN_in = in->numNodes;
117      out->numElements = 0;  #pragma omp parallel for
118      out->Id = NULL;      for (index_t e = 0; e < numElements; e++) {
119      out->Nodes = NULL;          const index_t k = index[e];
120      out->Tag = NULL;          Id[e] = in->Id[k];
121      out->Color = NULL;          Tag[e] = in->Tag[k];
122      out->minColor = 0;          Owner[e] = in->Owner[k];
123      out->maxColor = -1;          Color[e] = in->Color[k] + maxColor + 1;
124      out->jacobians = NULL;          for (int j = 0; j < std::min(numNodes, NN_in); j++)
125      out->jacobians_reducedQ = NULL;              Nodes[INDEX2(j, e, numNodes)] = in->Nodes[INDEX2(j, k, NN_in)];
   
     out->Owner = NULL;  
     out->numTagsInUse = 0;  
     out->tagsInUse = NULL;  
   
     out->MPIInfo = MPIInfo;  
   
     out->jacobians = Dudley_ElementFile_Jacobians_alloc();  
     out->jacobians_reducedQ = Dudley_ElementFile_Jacobians_alloc();  
   
     out->etype = etype;  
     out->numDim = Dims[out->etype];  
     out->numNodes = out->numDim + 1;  
     out->numLocalDim = localDims[out->etype];  
     out->numShapes = out->numLocalDim + 1;  
     out->ename = getElementName(out->etype);  
     return out;  
 }  
   
 /*  deallocates an element file: */  
 void Dudley_ElementFile_free(Dudley_ElementFile* in)  
 {  
     if (in != NULL)  
     {  
         Dudley_ElementFile_freeTable(in);  
         Dudley_ElementFile_Jacobians_dealloc(in->jacobians);  
         Dudley_ElementFile_Jacobians_dealloc(in->jacobians_reducedQ);  
         delete in;  
126      }      }
127        minColor = std::min(minColor, in->minColor + maxColor + 1);
128        maxColor = std::max(maxColor, in->maxColor + maxColor + 1);
129  }  }
130    
131  void Dudley_ElementFile_setElementDistribution(Dudley_ElementFile * in, dim_t * distribution)  void ElementFile::swapTable(ElementFile* other)
132  {  {
133      dim_t local_num_elements, e, num_elements = 0;      std::swap(numElements, other->numElements);
134      int myRank;      std::swap(Owner, other->Owner);
135      if (in == NULL) {      std::swap(Id, other->Id);
136          distribution[0] = num_elements;      std::swap(Nodes, other->Nodes);
137        std::swap(Tag, other->Tag);
138        std::swap(Color, other->Color);
139        std::swap(minColor, other->minColor);
140        std::swap(maxColor, other->maxColor);
141        std::swap(tagsInUse, other->tagsInUse);
142    }
143    
144    void ElementFile::optimizeOrdering()
145    {
146        if (numElements < 1)
147            return;
148    
149        util::ValueAndIndexList item_list(numElements);
150        index_t* index = new index_t[numElements];
151        ElementFile* out = new ElementFile(etype, MPIInfo);
152        out->allocTable(numElements);
153    
154    #pragma omp parallel for
155        for (index_t e = 0; e < numElements; e++) {
156            std::pair<index_t,index_t> entry(Nodes[INDEX2(0, e, numNodes)], e);
157            for (int i = 1; i < numNodes; i++)
158                entry.first = std::min(entry.first, Nodes[INDEX2(i, e, numNodes)]);
159            item_list[e] = entry;
160        }
161        util::sortValueAndIndex(item_list);
162    
163    #pragma omp parallel for
164        for (index_t e = 0; e < numElements; e++)
165            index[e] = item_list[e].second;
166    
167        out->gather(index, this);
168        swapTable(out);
169        delete out;
170        delete[] index;
171    }
172    
173    void ElementFile::setTags(int newTag, const escript::Data& mask)
174    {
175        const int numQuad = hasReducedIntegrationOrder(mask) ? 1 : numNodes;
176    
177        if (1 != mask.getDataPointSize()) {
178            throw DudleyException("ElementFile::setTags: number of components of mask must be 1.");
179        } else if (!mask.numSamplesEqual(numQuad, numElements)) {
180            throw DudleyException("ElementFile::setTags: illegal number of samples of mask Data object");
181        }
182    
183        if (mask.actsExpanded()) {
184    #pragma omp parallel for
185            for (index_t n = 0; n < numElements; n++) {
186                if (mask.getSampleDataRO(n)[0] > 0)
187                    Tag[n] = newTag;
188            }
189      } else {      } else {
190          if (in->MPIInfo->size > 1) {  #pragma omp parallel for
191              num_elements = 0;          for (index_t n = 0; n < numElements; n++) {
192              myRank = in->MPIInfo->rank;              const double* mask_array = mask.getSampleDataRO(n);
193  #pragma omp parallel private(local_num_elements)              bool check = false;
194              {              for (int q = 0; q < numQuad; q++)
195                  local_num_elements = 0;                  check = check || mask_array[q];
196  #pragma omp for private(e)              if (check)
197                  for (e = 0; e < in->numElements; e++) {                  Tag[n] = newTag;
                     if (in->Owner[e] == myRank)  
                         local_num_elements++;  
                 }  
 #pragma omp critical  
                 num_elements += local_num_elements;  
             }  
 #ifdef ESYS_MPI  
             MPI_Allgather(&num_elements, 1, MPI_INT, distribution, 1, MPI_INT, in->MPIInfo->comm);  
 #else  
             distribution[0] = num_elements;  
 #endif  
         } else {  
             distribution[0] = in->numElements;  
198          }          }
199      }      }
200        updateTagList();
201  }  }
202    
203  dim_t Dudley_ElementFile_getGlobalNumElements(Dudley_ElementFile* in)  void ElementFile::markNodes(std::vector<short>& mask, index_t offset) const
204  {  {
205      if (in == NULL)  #pragma omp parallel for
206          return 0;      for (index_t e = 0; e < numElements; e++) {
207            for (int i = 0; i < numNodes; i++) {
208                mask[Nodes[INDEX2(i, e, numNodes)] - offset] = 1;
209            }
210        }
211    }
212    
213      dim_t size, *distribution = NULL, out, p;  void ElementFile::relabelNodes(const index_t* newNode, index_t offset)
214      size = in->MPIInfo->size;  {
215      distribution = new  dim_t[size];  #pragma omp parallel for
216      Dudley_ElementFile_setElementDistribution(in, distribution);      for (index_t j = 0; j < numElements; j++) {
217      out = 0;          for (int i = 0; i < numNodes; i++) {
218      for (p = 0; p < size; ++p)              Nodes[INDEX2(i, j, numNodes)] =
219          out += distribution[p];                                newNode[Nodes[INDEX2(i, j, numNodes)] - offset];
220      delete[] distribution;          }
221      return out;      }
 }  
   
 dim_t Dudley_ElementFile_getMyNumElements(Dudley_ElementFile* in)  
 {  
     if (in == NULL)  
         return 0;  
   
     dim_t size, *distribution = NULL, out;  
     size = in->MPIInfo->size;  
     distribution = new  dim_t[size];  
     Dudley_ElementFile_setElementDistribution(in, distribution);  
     out = distribution[in->MPIInfo->rank];  
     delete[] distribution;  
     return out;  
 }  
   
 index_t Dudley_ElementFile_getFirstElement(Dudley_ElementFile* in)  
 {  
     if (in == NULL)  
         return 0;  
   
     dim_t size, *distribution = NULL, out, p;  
     size = in->MPIInfo->size;  
     distribution = new  dim_t[size];  
     Dudley_ElementFile_setElementDistribution(in, distribution);  
     out = 0;  
     for (p = 0; p < in->MPIInfo->rank; ++p)  
         out += distribution[p];  
     delete[] distribution;  
     return out;  
222  }  }
223    
224  } // namespace dudley  } // namespace dudley

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

  ViewVC Help
Powered by ViewVC 1.1.26