Actual source code: Completion.hh

  1: #ifndef included_ALE_Completion_hh
  2: #define included_ALE_Completion_hh

  4: #ifndef  included_ALE_Sections_hh
  5: #include <Sections.hh>
  6: #endif

  8: #ifndef  included_ALE_ParallelMapping_hh
  9: #include <ParallelMapping.hh>
 10: #endif

 12: #include <iostream>
 13: #include <fstream>

 15: namespace ALE {
 16:   class Completion {
 17:   public:
 18:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
 19:     static void completeSection(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
 20:       typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
 21:       //typedef ALE::Section<typename SendSection::point_type, typename SendSection::value_type> OverlapSection;
 22:       Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());

 24:       if (sendSection->debug()) {sendSection->view("Send Section");}
 25:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
 26:       if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
 27:       ALE::Pullback::InsertionBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
 28:       if (recvSection->debug()) {recvSection->view("Receieve Section");}
 29:     };
 30:     template<typename SendOverlap, typename RecvOverlap, typename SendSection, typename RecvSection>
 31:     static void completeSectionAdd(const Obj<SendOverlap>& sendOverlap, const Obj<RecvOverlap>& recvOverlap, const Obj<SendSection>& sendSection, const Obj<RecvSection>& recvSection) {
 32:       typedef ALE::Section<ALE::Pair<int, typename SendOverlap::source_type>, typename SendSection::value_type> OverlapSection;
 33:       Obj<OverlapSection> overlapSection = new OverlapSection(sendSection->comm(), sendSection->debug());

 35:       if (sendSection->debug()) {sendSection->view("Send Section");}
 36:       ALE::Pullback::SimpleCopy::copy(sendOverlap, recvOverlap, sendSection, overlapSection);
 37:       if (overlapSection->debug()) {overlapSection->view("Overlap Section");}
 38:       ALE::Pullback::AdditiveBinaryFusion::fuse(overlapSection, recvOverlap, recvSection);
 39:       if (recvSection->debug()) {recvSection->view("Receieve Section");}
 40:     };
 41:   };
 42:   namespace New {
 43:     template<typename Bundle_, typename Value_, typename Alloc_ = malloc_allocator<typename Bundle_::point_type> >
 44:     class Completion {
 45:     public:
 46:       typedef int                                                                         point_type;
 47:       typedef Value_                                                                      value_type;
 48:       typedef Bundle_                                                                     bundle_type;
 49:       typedef Alloc_                                                                      alloc_type;
 50:       typedef typename alloc_type::template rebind<int>::other                            int_alloc_type;
 51:       typedef typename alloc_type::template rebind<value_type>::other                     value_alloc_type;
 52:       typedef typename bundle_type::sieve_type                                            sieve_type;
 53:       typedef typename ALE::DiscreteSieve<point_type, alloc_type>                         dsieve_type;
 54:       typedef typename ALE::Topology<int, dsieve_type, alloc_type>                        topology_type;
 55:       typedef typename bundle_type::rank_type                                             rank_type;
 56:       typedef typename ALE::Sifter<point_type, rank_type, point_type>                     send_overlap_type;
 57:       typedef typename ALE::Sifter<rank_type, point_type, point_type>                     recv_overlap_type;
 58:       typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > send_sizer_type;
 59:       typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, int, int_alloc_type> > recv_sizer_type;
 60:       typedef typename ALE::New::ConeSizeSection<bundle_type, sieve_type>                 cone_size_section;
 61:       typedef typename ALE::New::ConeSection<sieve_type>                                  cone_section;
 62:       typedef typename ALE::New::SectionCompletion<bundle_type, value_type, alloc_type>   completion;
 63:     public:
 64:       template<typename PartitionType>
 65:       static void scatterSieve(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, const int dim, const Obj<sieve_type>& sieveNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const int height, const int numCells, const PartitionType assignment[]) {
 66:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
 67:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
 68:         int rank  = sieve->commRank();
 69:         int debug = sieve->debug();

 71:         // Create local sieve
 72:         const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(height);
 73:         int e = 0;

 75:         if (sieve->debug()) {
 76:           int e2 = 0;
 77:           for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
 78:             std::cout << "assignment["<<*e_iter<<"]" << assignment[e2++] << std::endl;
 79:           }
 80:         }
 81:         PetscTruth flg;
 82:         PetscOptionsHasName(PETSC_NULL, "-output_partition", &flg);
 83:         if (flg) {
 84:           ostringstream fname;
 85:           fname << "part." << sieve->commSize() << ".dat";
 86:           std::ofstream f(fname.str().c_str());
 87:           int e2 = 0;
 88:           f << sieve->commSize() << std::endl;
 89:           for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
 90:             f << assignment[e2++] << std::endl;
 91:           }
 92:         }
 93:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
 94:           if (assignment[e] == rank) {
 95:             Obj<typename sieve_type::coneSet> current = new typename sieve_type::coneSet();
 96:             Obj<typename sieve_type::coneSet> next    = new typename sieve_type::coneSet();
 97:             Obj<typename sieve_type::coneSet> tmp;

 99:             current->insert(*e_iter);
100:             while(current->size()) {
101:               for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
102:                 const Obj<typename sieve_type::traits::coneSequence>& cone = sieve->cone(*p_iter);
103: 
104:                 for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
105:                   sieveNew->addArrow(*c_iter, *p_iter, c_iter.color());
106:                   next->insert(*c_iter);
107:                 }
108:               }
109:               tmp = current; current = next; next = tmp;
110:               next->clear();
111:             }
112:             if (height) {
113:               current->insert(*e_iter);
114:               while(current->size()) {
115:                 for(typename sieve_type::coneSet::const_iterator p_iter = current->begin(); p_iter != current->end(); ++p_iter) {
116:                   const Obj<typename sieve_type::traits::supportSequence>& support = sieve->support(*p_iter);
117: 
118:                   for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
119:                     sieveNew->addArrow(*p_iter, *s_iter, s_iter.color());
120:                     next->insert(*s_iter);
121:                   }
122:                 }
123:                 tmp = current; current = next; next = tmp;
124:                 next->clear();
125:               }
126:             }
127:           }
128:           e++;
129:         }
130:         // Complete sizer section
131:         typedef typename ALE::New::PartitionSizeSection<bundle_type, PartitionType> partition_size_section;
132:         typedef typename ALE::New::PartitionSection<bundle_type, PartitionType>     partition_section;
133:         Obj<topology_type>          secTopology          = completion::createSendTopology(sendOverlap);
134:         Obj<partition_size_section> partitionSizeSection = new partition_size_section(bundle, height, numCells, assignment);
135:         Obj<partition_section>      partitionSection     = new partition_section(bundle, height, numCells, assignment);
136:         Obj<send_section_type>      sendSection          = new send_section_type(sieve->comm(), sieve->debug());
137:         Obj<recv_section_type>      recvSection          = new recv_section_type(sieve->comm(), sendSection->getTag(), sieve->debug());

139:         completion::completeSection(sendOverlap, recvOverlap, partitionSizeSection, partitionSection, sendSection, recvSection);
140:         // Unpack the section into the overlap
141:         sendOverlap->clear();
142:         recvOverlap->clear();
143:         const typename send_section_type::sheaf_type& sendPatches = sendSection->getPatches();

145:         for(typename send_section_type::sheaf_type::const_iterator p_iter = sendPatches.begin(); p_iter != sendPatches.end(); ++p_iter) {
146:           const typename send_section_type::patch_type               rank    = p_iter->first;
147:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
148:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

150:           for(typename send_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
151:             const typename send_section_type::value_type *points = section->restrictPoint(*c_iter);
152:             int                                           size   = section->getFiberDimension(*c_iter);

154:             for(int p = 0; p < size; p++) {
155:               sendOverlap->addArrow(points[p], rank, points[p]);
156:             }
157:           }
158:         }
159:         const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();

161:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
162:           const typename send_section_type::patch_type               rank    = p_iter->first;
163:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
164:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

166:           for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
167:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
168:             int                                           size   = section->getFiberDimension(*c_iter);

170:             for(int p = 0; p < size; p++) {
171:               recvOverlap->addArrow(rank, points[p], points[p]);
172:             }
173:           }
174:         }
175:         if (debug) {
176:           sendOverlap->view(std::cout, "Send overlap for points");
177:           recvOverlap->view(std::cout, "Receive overlap for points");
178:         }
179:         // Receive the point section
180:         ALE::New::Completion<bundle_type, value_type>::scatterCones(sieve, sieveNew, sendOverlap, recvOverlap, bundle, height);
181:         if (height) {
182:           ALE::New::Completion<bundle_type, value_type>::scatterSupports(sieve, sieveNew, sendOverlap, recvOverlap, bundle, bundle->depth()-height);
183:         }
184:       };
185:       template<typename SifterType>
186:       static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumHeight = 0) {
187:         typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
188:         typedef typename ALE::New::ConeSection<SifterType>                  cone_section;
189:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
190:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
191:         Obj<topology_type>     secTopology     = completion::createSendTopology(sendOverlap);
192:         Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter, minimumHeight);
193:         Obj<cone_section>      coneSection     = new cone_section(sifter);
194:         Obj<send_section_type> sendSection     = new send_section_type(sifter->comm(), sifter->debug());
195:         Obj<recv_section_type> recvSection     = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

197:         completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
198:         // Unpack the section into the sieve
199:         const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();

201:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
202:           const Obj<typename recv_section_type::section_type>&        section = p_iter->second;
203:           const typename recv_section_type::section_type::chart_type& chart   = section->getChart();

205:           for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
206:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
207:             int size = section->getFiberDimension(*c_iter);
208:             int c    = 0;

210:             for(int p = 0; p < size; p++) {
211:               sifterNew->addArrow(points[p], *c_iter, c++);
212:             }
213:           }
214:         }
215:       };
216:       template<typename SifterType, typename Renumbering>
217:       static void scatterCones(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, Renumbering& renumbering, const Obj<bundle_type>& bundle = NULL) {
218:         typedef typename ALE::New::ConeSizeSection<bundle_type, SifterType> cone_size_section;
219:         typedef typename ALE::New::ConeSection<SifterType>                  cone_section;
220:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
221:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
222:         Obj<topology_type>     secTopology     = completion::createSendTopology(sendOverlap);
223:         Obj<cone_size_section> coneSizeSection = new cone_size_section(bundle, sifter);
224:         Obj<cone_section>      coneSection     = new cone_section(sifter);
225:         Obj<send_section_type> sendSection     = new send_section_type(sifter->comm(), sifter->debug());
226:         Obj<recv_section_type> recvSection     = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

228:         completion::completeSection(sendOverlap, recvOverlap, coneSizeSection, coneSection, sendSection, recvSection);
229:         // Unpack the section into the sieve
230:         const typename recv_section_type::sheaf_type& patches = recvSection->getPatches();

232:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = patches.begin(); p_iter != patches.end(); ++p_iter) {
233:           const Obj<typename recv_section_type::section_type>&        section = p_iter->second;
234:           const typename recv_section_type::section_type::chart_type& chart   = section->getChart();

236:           for(typename recv_section_type::section_type::chart_type::const_iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
237:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
238:             int size = section->getFiberDimension(*c_iter);
239:             int c    = 0;

241:             for(int p = 0; p < size; p++) {
242:               sifterNew->addArrow(points[p], renumbering[*c_iter], c++);
243:             }
244:           }
245:         }
246:       };
247:       template<typename SifterType>
248:       static void scatterSupports(const Obj<SifterType>& sifter, const Obj<SifterType>& sifterNew, const Obj<send_overlap_type>& sendOverlap, const Obj<recv_overlap_type>& recvOverlap, const Obj<bundle_type>& bundle = NULL, const int minimumDepth = 0) {
249:         typedef typename ALE::New::SupportSizeSection<bundle_type, SifterType> support_size_section;
250:         typedef typename ALE::New::SupportSection<SifterType>                  support_section;
251:         typedef typename ALE::Field<send_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > send_section_type;
252:         typedef typename ALE::Field<recv_overlap_type, int, ALE::Section<point_type, value_type, value_alloc_type> > recv_section_type;
253:         Obj<topology_type>        secTopology        = completion::createSendTopology(sendOverlap);
254:         Obj<support_size_section> supportSizeSection = new support_size_section(bundle, sifter, minimumDepth);
255:         Obj<support_section>      supportSection     = new support_section(sifter);
256:         Obj<send_section_type>    sendSection        = new send_section_type(sifter->comm(), sifter->debug());
257:         Obj<recv_section_type>    recvSection        = new recv_section_type(sifter->comm(), sendSection->getTag(), sifter->debug());

259:         completion::completeSection(sendOverlap, recvOverlap, supportSizeSection, supportSection, sendSection, recvSection);
260:         // Unpack the section into the sieve
261:         const typename recv_section_type::sheaf_type& recvPatches = recvSection->getPatches();

263:         for(typename recv_section_type::sheaf_type::const_iterator p_iter = recvPatches.begin(); p_iter != recvPatches.end(); ++p_iter) {
264:           const Obj<typename send_section_type::section_type>&       section = p_iter->second;
265:           const typename send_section_type::section_type::chart_type chart   = section->getChart();

267:           for(typename recv_section_type::section_type::chart_type::iterator c_iter = chart.begin(); c_iter != chart.end(); ++c_iter) {
268:             const typename recv_section_type::value_type *points = section->restrictPoint(*c_iter);
269:             int                                           size   = section->getFiberDimension(*c_iter);
270:             int                                           c      = 0;

272:             for(int p = 0; p < size; p++) {
273:               sifterNew->addArrow(*c_iter, points[p], c++);
274:             }
275:           }
276:         }
277:       };
278:     };

280:     template<typename Value_>
281:     class ParallelFactory {
282:     public:
283:       typedef Value_ value_type;
284:     protected:
285:       int          _debug;
286:       MPI_Datatype _mpiType;
287:     protected:
288:       MPI_Datatype constructMPIType() {
289:         if (sizeof(value_type) == 1) {
290:           return MPI_BYTE;
291:         } else if (sizeof(value_type) == 2) {
292:           return MPI_SHORT;
293:         } else if (sizeof(value_type) == 4) {
294:           return MPI_INT;
295:         } else if (sizeof(value_type) == 8) {
296:           return MPI_DOUBLE;
297:         } else if (sizeof(value_type) == 28) {
298:           int          blen[2];
299:           MPI_Aint     indices[2];
300:           MPI_Datatype oldtypes[2], newtype;
301:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_INT;
302:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
303:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
304:           MPI_Type_commit(&newtype);
305:           return newtype;
306:         } else if (sizeof(value_type) == 32) {
307:           int          blen[2];
308:           MPI_Aint     indices[2];
309:           MPI_Datatype oldtypes[2], newtype;
310:           blen[0] = 1; indices[0] = 0;           oldtypes[0] = MPI_DOUBLE;
311:           blen[1] = 3; indices[1] = sizeof(int); oldtypes[1] = MPI_DOUBLE;
312:           MPI_Type_struct(2, blen, indices, oldtypes, &newtype);
313:           MPI_Type_commit(&newtype);
314:           return newtype;
315:         }
316:         throw ALE::Exception("Cannot determine MPI type for value type");
317:       };
318:       ParallelFactory(const int debug) : _debug(debug) {
319:         this->_mpiType = this->constructMPIType();
320:       };
321:     public:
322:       ~ParallelFactory() {};
323:     public:
324:       static const Obj<ParallelFactory>& singleton(const int debug, bool cleanup = false) {
325:         static Obj<ParallelFactory> *_singleton = NULL;

327:         if (cleanup) {
328:           if (debug) {std::cout << "Destroying ParallelFactory" << std::endl;}
329:           if (_singleton) {delete _singleton;}
330:           _singleton = NULL;
331:         } else if (_singleton == NULL) {
332:           if (debug) {std::cout << "Creating new ParallelFactory" << std::endl;}
333:           _singleton  = new Obj<ParallelFactory>();
334:           *_singleton = new ParallelFactory(debug);
335:         }
336:         return *_singleton;
337:       };
338:     public: // Accessors
339:       int debug() const {return this->_debug;};
340:       MPI_Datatype getMPIType() const {return this->_mpiType;};
341:     };
342:   }
343: }
344: #endif