Actual source code: Sections.hh

  1: #ifndef included_ALE_Sections_hh
  2: #define included_ALE_Sections_hh

  4: #ifndef  included_ALE_Numbering_hh
  5: #include <Numbering.hh>
  6: #endif

  8: namespace ALE {
  9:   template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
 10:   class BaseSection : public ALE::ParallelObject {
 11:   public:
 12:     typedef Sieve_                                    sieve_type;
 13:     typedef Alloc_                                    alloc_type;
 14:     typedef int                                       value_type;
 15:     typedef typename sieve_type::target_type          point_type;
 16:     typedef typename sieve_type::traits::baseSequence chart_type;
 17:   protected:
 18:     Obj<sieve_type> _sieve;
 19:     chart_type      _chart;
 20:     int             _sizes[2];
 21:   public:
 22:     BaseSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
 23:     ~BaseSection() {};
 24:   public: // Verifiers
 25:     bool hasPoint(const point_type& point) const {
 26:       return this->_sieve->baseContains(point);
 27:     };
 28:   public:
 29:     const chart_type& getChart() const {
 30:       return this->_chart;
 31:     };
 32:     const int getFiberDimension(const point_type& p) const {
 33:       return this->hasPoint(p) ? 1 : 0;
 34:     };
 35:     const value_type *restrictSpace() const {
 36:       return this->_sizes;
 37:     };
 38:     const value_type *restrictPoint(const point_type& p) const {
 39:       if (this->hasPoint(p)) return this->_sizes;
 40:       return &this->_sizes[1];
 41:     };
 42:   };

 44:   template<typename Sieve_, typename Alloc_ = malloc_allocator<int> >
 45:   class ConeSizeSection : public ALE::ParallelObject {
 46:   public:
 47:     typedef Sieve_                              sieve_type;
 48:     typedef Alloc_                              alloc_type;
 49:     typedef int                                 value_type;
 50:     typedef typename sieve_type::target_type    point_type;
 51:     typedef BaseSection<sieve_type, alloc_type> atlas_type;
 52:     typedef typename atlas_type::chart_type     chart_type;
 53:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
 54:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
 55:   protected:
 56:     Obj<sieve_type> _sieve;
 57:     Obj<atlas_type> _atlas;
 58:     int             _size;
 59:   public:
 60:     ConeSizeSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
 61:       atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
 62:       atlas_alloc_type().construct(pAtlas, atlas_type(sieve));
 63:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
 64:     };
 65:     ~ConeSizeSection() {};
 66:   public: // Verifiers
 67:     bool hasPoint(const point_type& point) {
 68:       return this->_atlas->hasPoint(point);
 69:     };
 70:   public: // Accessors
 71:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
 72:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
 73:   public:
 74:     const int getFiberDimension(const point_type& p) {
 75:       return this->hasPoint(p) ? 1 : 0;
 76:     };
 77:     const value_type *restrictPoint(const point_type& p) {
 78:       this->_size = this->_sieve->cone(p)->size();
 79:       return &this->_size;
 80:     };
 81:   };

 83:   template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::source_type> >
 84:   class ConeSection : public ALE::ParallelObject {
 85:   public:
 86:     typedef Sieve_                                  sieve_type;
 87:     typedef Alloc_                                  alloc_type;
 88:     typedef typename sieve_type::target_type        point_type;
 89:     typedef typename sieve_type::source_type        value_type;
 90:     typedef ConeSizeSection<sieve_type, alloc_type> atlas_type;
 91:     typedef typename atlas_type::chart_type         chart_type;
 92:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
 93:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
 94:   protected:
 95:     Obj<sieve_type> _sieve;
 96:     Obj<atlas_type> _atlas;
 97:     alloc_type      _allocator;
 98:   public:
 99:     ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
100:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
101:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
102:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
103:     };
104:     ~ConeSection() {};
105:   protected:
106:     value_type *getRawArray(const int size) {
107:       static value_type *array   = NULL;
108:       static int         maxSize = 0;

110:       if (size > maxSize) {
111:         const value_type dummy(0);

113:         if (array) {
114:           for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
115:           this->_allocator.deallocate(array, maxSize);
116:         }
117:         maxSize = size;
118:         array   = this->_allocator.allocate(maxSize);
119:         for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
120:       };
121:       return array;
122:     };
123:   public: // Verifiers
124:     bool hasPoint(const point_type& point) {
125:       return this->_atlas->hasPoint(point);
126:     };
127:   public: // Accessors
128:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
129:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
130:   public: // Sizes and storage
131:     int getFiberDimension(const point_type& p) {
132:       return this->_atlas->restrictPoint(p)[0];
133:     };
134:   public: // Restriction and update
135:     const value_type *restrictPoint(const point_type& p) {
136:       const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
137:       value_type *array = this->getRawArray(cone->size());
138:       int         c     = 0;

140:       for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
141:         array[c++] = *c_iter;
142:       }
143:       return array;
144:     };
145:   };

147:   template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
148:   class BaseSectionV : public ALE::ParallelObject {
149:   public:
150:     typedef Sieve_                                    sieve_type;
151:     typedef Alloc_                                    alloc_type;
152:     typedef int                                       value_type;
153:     typedef typename sieve_type::target_type          point_type;
154:     //typedef typename sieve_type::traits::baseSequence chart_type;
155:     typedef int chart_type;
156:   protected:
157:     Obj<sieve_type> _sieve;
158:     //chart_type      _chart;
159:     int             _sizes[2];
160:   public:
161:     //BaseSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
162:     BaseSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {_sizes[0] = 1; _sizes[1] = 0;};
163:     ~BaseSectionV() {};
164:   public: // Verifiers
165:     bool hasPoint(const point_type& point) const {
166:       return this->_sieve->baseContains(point);
167:     };
168:   public:
169:     //const chart_type& getChart() const {
170:     //  return this->_chart;
171:     //};
172:     const int getFiberDimension(const point_type& p) const {
173:       return this->hasPoint(p) ? 1 : 0;
174:     };
175:     const value_type *restrictSpace() const {
176:       return this->_sizes;
177:     };
178:     const value_type *restrictPoint(const point_type& p) const {
179:       if (this->hasPoint(p)) return this->_sizes;
180:       return &this->_sizes[1];
181:     };
182:   };

184:   template<typename Sieve_, typename Alloc_ = malloc_allocator<int> >
185:   class ConeSizeSectionV : public ALE::ParallelObject {
186:   public:
187:     typedef Sieve_                               sieve_type;
188:     typedef Alloc_                               alloc_type;
189:     typedef int                                  value_type;
190:     typedef typename sieve_type::target_type     point_type;
191:     typedef BaseSectionV<sieve_type, alloc_type> atlas_type;
192:     typedef typename atlas_type::chart_type      chart_type;
193:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
194:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
195:   protected:
196:     Obj<sieve_type> _sieve;
197:     Obj<atlas_type> _atlas;
198:     int             _size;
199:   public:
200:     ConeSizeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
201:       atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
202:       atlas_alloc_type().construct(pAtlas, atlas_type(sieve));
203:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
204:     };
205:     ~ConeSizeSectionV() {};
206:   public: // Verifiers
207:     bool hasPoint(const point_type& point) {
208:       return this->_atlas->hasPoint(point);
209:     };
210:   public: // Accessors
211:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
212:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
213:   public:
214:     const int getFiberDimension(const point_type& p) {
215:       return this->hasPoint(p) ? 1 : 0;
216:     };
217:     const value_type *restrictPoint(const point_type& p) {
218:       this->_size = this->_sieve->getConeSize(p);
219:       return &this->_size;
220:     };
221:   };

223:   template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::source_type> >
224:   class ConeSectionV : public ALE::ParallelObject {
225:   public:
226:     typedef Sieve_                                   sieve_type;
227:     typedef Alloc_                                   alloc_type;
228:     typedef typename sieve_type::target_type         point_type;
229:     typedef typename sieve_type::source_type         value_type;
230:     typedef ConeSizeSectionV<sieve_type, alloc_type> atlas_type;
231:     typedef typename atlas_type::chart_type          chart_type;
232:     typedef typename ISieveVisitor::PointRetriever<sieve_type> visitor_type;
233:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
234:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
235:   protected:
236:     Obj<sieve_type> _sieve;
237:     Obj<atlas_type> _atlas;
238:     visitor_type   *_cV;
239:     alloc_type      _allocator;
240:   public:
241:     ConeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
242:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
243:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
244:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
245:       this->_cV        = new visitor_type(std::max(0, sieve->getMaxConeSize()));
246:     };
247:     ~ConeSectionV() {
248:       delete this->_cV;
249:     };
250:   public: // Verifiers
251:     bool hasPoint(const point_type& point) {
252:       return this->_atlas->hasPoint(point);
253:     };
254:   public: // Accessors
255:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
256:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
257:   public: // Sizes and storage
258:     int getFiberDimension(const point_type& p) {
259:       return this->_atlas->restrictPoint(p)[0];
260:     };
261:   public: // Restriction and update
262:     const value_type *restrictPoint(const point_type& p) {
263:       this->_cV->clear();
264:       this->_sieve->cone(p, *this->_cV);
265:       return this->_cV->getPoints();
266:     };
267:   };

269:   template<typename Sieve_, typename Alloc_ = malloc_allocator<OrientedPoint<typename Sieve_::source_type> > >
270:   class OrientedConeSectionV : public ALE::ParallelObject {
271:   public:
272:     typedef Sieve_                                   sieve_type;
273:     typedef Alloc_                                   alloc_type;
274:     typedef typename sieve_type::target_type         point_type;
275:     typedef OrientedPoint<typename sieve_type::source_type> value_type;
276:     typedef typename alloc_type::template rebind<int>::other int_alloc_type;
277:     typedef ConeSizeSectionV<sieve_type, int_alloc_type> atlas_type;
278:     typedef typename atlas_type::chart_type          chart_type;
279:     typedef typename ISieveVisitor::PointRetriever<sieve_type> visitor_type;
280:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
281:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
282:   protected:
283:     Obj<sieve_type> _sieve;
284:     Obj<atlas_type> _atlas;
285:     visitor_type   *_cV;
286:     alloc_type      _allocator;
287:   public:
288:     OrientedConeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
289:       atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
290:       atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
291:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
292:       this->_cV        = new visitor_type(std::max(0, sieve->getMaxConeSize()));
293:     };
294:     ~OrientedConeSectionV() {
295:       delete this->_cV;
296:     };
297:   public: // Verifiers
298:     bool hasPoint(const point_type& point) {
299:       return this->_atlas->hasPoint(point);
300:     };
301:   public: // Accessors
302:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
303:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
304:   public: // Sizes and storage
305:     int getFiberDimension(const point_type& p) {
306:       return this->_atlas->restrictPoint(p)[0];
307:     };
308:   public: // Restriction and update
309:     const value_type *restrictPoint(const point_type& p) {
310:       this->_cV->clear();
311:       this->_sieve->orientedCone(p, *this->_cV);
312:       return (const value_type *) this->_cV->getOrientedPoints();
313:     };
314:   };

316:   template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
317:   class LabelBaseSection : public ALE::ParallelObject {
318:   public:
319:     typedef Sieve_                                    sieve_type;
320:     typedef Label_                                    label_type;
321:     typedef Alloc_                                    alloc_type;
322:     typedef int                                       value_type;
323:     typedef typename sieve_type::target_type          point_type;
324:     typedef typename sieve_type::traits::baseSequence chart_type;
325:   protected:
326:     Obj<sieve_type> _sieve;
327:     Obj<label_type> _label;
328:     chart_type      _chart;
329:     int             _sizes[2];
330:   public:
331:     LabelBaseSection(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
332:     ~LabelBaseSection() {};
333:   public: // Verifiers
334:     bool hasPoint(const point_type& point) const {
335:       return this->_label->cone(point)->size() ? true : false;
336:     };
337:   public:
338:     const chart_type& getChart() const {
339:       return this->_chart;
340:     };
341:     const int getFiberDimension(const point_type& p) const {
342:       return this->hasPoint(p) ? 1 : 0;
343:     };
344:     const value_type *restrictSpace() const {
345:       return this->_sizes;
346:     };
347:     const value_type *restrictPoint(const point_type& p) const {
348:       if (this->hasPoint(p)) return this->_sizes;
349:       return &this->_sizes[1];
350:     };
351:   };

353:   template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<int> >
354:   class LabelSection : public ALE::ParallelObject {
355:   public:
356:     typedef Sieve_                              sieve_type;
357:     typedef Label_                              label_type;
358:     typedef Alloc_                              alloc_type;
359:     typedef int                                 value_type;
360:     typedef typename sieve_type::target_type    point_type;
361:     typedef LabelBaseSection<sieve_type, label_type, alloc_type> atlas_type;
362:     typedef typename atlas_type::chart_type     chart_type;
363:     typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
364:     typedef typename atlas_alloc_type::pointer                      atlas_ptr;
365:   protected:
366:     Obj<sieve_type> _sieve;
367:     Obj<label_type> _label;
368:     Obj<atlas_type> _atlas;
369:     int             _size;
370:     int             _value;
371:   public:
372:     LabelSection(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label) {
373:       atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
374:       atlas_alloc_type().construct(pAtlas, atlas_type(sieve, label));
375:       this->_atlas     = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
376:     };
377:     ~LabelSection() {};
378:   public: // Verifiers
379:     bool hasPoint(const point_type& point) {
380:       return this->_atlas->hasPoint(point);
381:     };
382:   public: // Accessors
383:     const Obj<atlas_type>& getAtlas() {return this->_atlas;};
384:     void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
385:   public:
386:     const int getFiberDimension(const point_type& p) {
387:       return this->hasPoint(p) ? 1 : 0;
388:     };
389:     const value_type *restrictPoint(const point_type& p) {
390:       this->_value = *this->_label->cone(p)->begin();
391:       return &this->_value;
392:     };
393:   };

395:   namespace New {
396:     // This section takes an existing section, and reports instead the fiber dimensions as values
397:     template<typename Section_>
398:     class SizeSection : public ALE::ParallelObject {
399:     public:
400:       typedef Section_                          section_type;
401:       typedef typename section_type::point_type point_type;
402:       typedef int                               value_type;
403:     protected:
404:       Obj<section_type> _section;
405:       value_type        _size;
406:     public:
407:       SizeSection(const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, section->debug()), _section(section) {};
408:       virtual ~SizeSection() {};
409:     public:
410:       bool hasPoint(const point_type& point) {
411:         return this->_section->hasPoint(point);
412:       };
413:       const value_type *restrictPoint(const point_type& p) {
414:         this->_size = this->_section->getFiberDimension(p);
415:         return &this->_size;
416:       };
417:     public:
418:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
419:         this->_section->view(name, comm);
420:       };
421:     };

423:     // This section reports as values the size of the partition associated with the partition point
424:     template<typename Bundle_, typename Marker_>
425:     class PartitionSizeSection : public ALE::ParallelObject {
426:     public:
427:       typedef Bundle_                          bundle_type;
428:       typedef typename bundle_type::sieve_type sieve_type;
429:       typedef typename bundle_type::point_type point_type;
430:       typedef Marker_                          marker_type;
431:       typedef int                              value_type;
432:       typedef std::map<marker_type, int>       sizes_type;
433:     protected:
434:       sizes_type _sizes;
435:       int        _height;
436:       void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
437:         const Obj<typename bundle_type::label_sequence>& cells      = bundle->heightStratum(this->_height);
438:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
439:         std::map<marker_type, std::set<point_type> >     points;

441:         if (numElements != (int) cells->size()) {
442:           throw ALE::Exception("Partition size does not match the number of elements");
443:         }
444:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
445:           typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
446:           const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
447:           const int idx = cNumbering->getIndex(*e_iter);

449:           points[partition[idx]].insert(closure->begin(), closure->end());
450:           if (this->_height > 0) {
451:             const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);

453:             points[partition[idx]].insert(star->begin(), star->end());
454:           }
455:         }
456:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
457:           this->_sizes[p_iter->first] = p_iter->second.size();
458:         }
459:       };
460:     public:
461:       PartitionSizeSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
462:         this->_init(bundle, numElements, partition);
463:       };
464:       virtual ~PartitionSizeSection() {};
465:     public:
466:       bool hasPoint(const point_type& point) {return true;};
467:       const value_type *restrictPoint(const point_type& p) {
468:         return &this->_sizes[p];
469:       };
470:     public:
471:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
472:         ostringstream txt;
473:         int rank;

475:         if (comm == MPI_COMM_NULL) {
476:           comm = this->comm();
477:           rank = this->commRank();
478:         } else {
479:           MPI_Comm_rank(comm, &rank);
480:         }
481:         if (name == "") {
482:           if(rank == 0) {
483:             txt << "viewing a PartitionSizeSection" << std::endl;
484:           }
485:         } else {
486:           if(rank == 0) {
487:             txt << "viewing PartitionSizeSection '" << name << "'" << std::endl;
488:           }
489:         }
490:         for(typename sizes_type::const_iterator s_iter = this->_sizes.begin(); s_iter != this->_sizes.end(); ++s_iter) {
491:           const marker_type& partition = s_iter->first;
492:           const value_type   size      = s_iter->second;

494:           txt << "[" << this->commRank() << "]: Partition " << partition << " size " << size << std::endl;
495:         }
496:         PetscSynchronizedPrintf(comm, txt.str().c_str());
497:         PetscSynchronizedFlush(comm);
498:       };
499:     };

501:     template<typename Point_>
502:     class PartitionDomain {
503:     public:
504:       typedef Point_ point_type;
505:     public:
506:       PartitionDomain() {};
507:       ~PartitionDomain() {};
508:     public:
509:       int count(const point_type& point) const {return 1;};
510:     };

512:     // This section returns the points in each partition
513:     template<typename Bundle_, typename Marker_>
514:     class PartitionSection : public ALE::ParallelObject {
515:     public:
516:       typedef Bundle_                            bundle_type;
517:       typedef typename bundle_type::sieve_type   sieve_type;
518:       typedef typename bundle_type::point_type   point_type;
519:       typedef Marker_                            marker_type;
520:       typedef int                                value_type;
521:       typedef std::map<marker_type, point_type*> points_type;
522:       typedef PartitionDomain<point_type>        chart_type;
523:     protected:
524:       points_type _points;
525:       chart_type  _domain;
526:       int         _height;
527:       void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
528:         // Should check for patch 0
529:         const Obj<typename bundle_type::label_sequence>& cells      = bundle->heightStratum(this->_height);
530:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
531:         std::map<marker_type, std::set<point_type> >     points;
532:         std::map<marker_type, int>                       offsets;

534:         if (numElements != (int) cells->size()) {
535:           throw ALE::Exception("Partition size does not match the number of elements");
536:         }
537:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
538:           typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
539:           const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
540:           const int idx = cNumbering->getIndex(*e_iter);

542:           points[partition[idx]].insert(closure->begin(), closure->end());
543:           if (this->_height > 0) {
544:             const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);

546:             points[partition[idx]].insert(star->begin(), star->end());
547:           }
548:         }
549:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
550:           this->_points[p_iter->first] = new point_type[p_iter->second.size()];
551:           offsets[p_iter->first] = 0;
552:           for(typename std::set<point_type>::const_iterator s_iter = p_iter->second.begin(); s_iter != p_iter->second.end(); ++s_iter) {
553:             this->_points[p_iter->first][offsets[p_iter->first]++] = *s_iter;
554:           }
555:         }
556:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
557:           if (offsets[p_iter->first] != (int) p_iter->second.size()) {
558:             ostringstream txt;
559:             txt << "Invalid offset for partition " << p_iter->first << ": " << offsets[p_iter->first] << " should be " << p_iter->second.size();
560:             throw ALE::Exception(txt.str().c_str());
561:           }
562:         }
563:       };
564:     public:
565:       PartitionSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
566:         this->_init(bundle, numElements, partition);
567:       };
568:       virtual ~PartitionSection() {
569:         for(typename points_type::iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
570:           delete [] p_iter->second;
571:         }
572:       };
573:     public:
574:       const chart_type& getChart() {return this->_domain;};
575:       bool hasPoint(const point_type& point) {return true;};
576:       const value_type *restrictPoint(const point_type& p) {
577:         return this->_points[p];
578:       };
579:     public:
580:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
581:         ostringstream txt;
582:         int rank;

584:         if (comm == MPI_COMM_NULL) {
585:           comm = this->comm();
586:           rank = this->commRank();
587:         } else {
588:           MPI_Comm_rank(comm, &rank);
589:         }
590:         if (name == "") {
591:           if(rank == 0) {
592:             txt << "viewing a PartitionSection" << std::endl;
593:           }
594:         } else {
595:           if(rank == 0) {
596:             txt << "viewing PartitionSection '" << name << "'" << std::endl;
597:           }
598:         }
599:         for(typename points_type::const_iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
600:           const marker_type& partition  = p_iter->first;
601:           //const point_type *points = p_iter->second;

603:           txt << "[" << this->commRank() << "]: Partition " << partition << std::endl;
604:         }
605:         if (this->_points.size() == 0) {
606:           txt << "[" << this->commRank() << "]: empty" << std::endl;
607:         }
608:         PetscSynchronizedPrintf(comm, txt.str().c_str());
609:         PetscSynchronizedFlush(comm);
610:       };
611:     };

613:     template<typename Bundle_, typename Sieve_>
614:     class ConeSizeSection : public ALE::ParallelObject {
615:     public:
616:       typedef ConeSizeSection<Bundle_, Sieve_> section_type;
617:       typedef int                              patch_type;
618:       typedef Bundle_                          bundle_type;
619:       typedef Sieve_                           sieve_type;
620:       typedef typename bundle_type::point_type point_type;
621:       typedef int                              value_type;
622:     protected:
623:       Obj<bundle_type> _bundle;
624:       Obj<sieve_type>  _sieve;
625:       value_type       _size;
626:       int              _minHeight;
627:       Obj<section_type> _section;
628:     public:
629:       ConeSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumHeight = 0) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _bundle(bundle), _sieve(sieve), _minHeight(minimumHeight) {
630:         this->_section = this;
631:         this->_section.addRef();
632:       };
633:       virtual ~ConeSizeSection() {};
634:     public: // Verifiers
635:       bool hasPoint(const point_type& point) {return true;};
636:     public: // Restriction
637:       const value_type *restrictPoint(const point_type& p) {
638:         if ((this->_minHeight == 0) || (this->_bundle->height(p) >= this->_minHeight)) {
639:           this->_size = this->_sieve->cone(p)->size();
640:         } else {
641:           this->_size = 0;
642:         }
643:         return &this->_size;
644:       };
645:     public: // Adapter
646:       const Obj<section_type>& getSection(const patch_type& patch) {
647:         return this->_section;
648:       };
649:     public:
650:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
651:         ostringstream txt;
652:         int rank;

654:         if (comm == MPI_COMM_NULL) {
655:           comm = this->comm();
656:           rank = this->commRank();
657:         } else {
658:           MPI_Comm_rank(comm, &rank);
659:         }
660:         if (name == "") {
661:           if(rank == 0) {
662:             txt << "viewing a ConeSizeSection" << std::endl;
663:           }
664:         } else {
665:           if(rank == 0) {
666:             txt << "viewing ConeSizeSection '" << name << "'" << std::endl;
667:           }
668:         }
669:         PetscSynchronizedPrintf(comm, txt.str().c_str());
670:         PetscSynchronizedFlush(comm);
671:       };
672:     };

674:     template<typename Sieve_>
675:     class ConeSection : public ALE::ParallelObject {
676:     public:
677:       typedef Sieve_                           sieve_type;
678:       typedef typename sieve_type::target_type point_type;
679:       typedef typename sieve_type::source_type value_type;
680:       typedef PartitionDomain<sieve_type>      chart_type;
681:     protected:
682:       Obj<sieve_type> _sieve;
683:       int             _coneSize;
684:       value_type     *_cone;
685:       chart_type      _domain;
686:       void ensureCone(const int size) {
687:         if (size > this->_coneSize) {
688:           if (this->_cone) delete [] this->_cone;
689:           this->_coneSize = size;
690:           this->_cone     = new value_type[this->_coneSize];
691:         }
692:       };
693:     public:
694:       ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _coneSize(-1), _cone(NULL) {};
695:       virtual ~ConeSection() {if (this->_cone) delete [] this->_cone;};
696:     public:
697:       const chart_type& getChart() {return this->_domain;};
698:       bool hasPoint(const point_type& point) {return true;};
699:       const value_type *restrictPoint(const point_type& p) {
700:         const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
701:         int c = 0;

703:         this->ensureCone(cone->size());
704:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
705:           this->_cone[c++] = *c_iter;
706:         }
707:         return this->_cone;
708:       };
709:     public:
710:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
711:         ostringstream txt;
712:         int rank;

714:         if (comm == MPI_COMM_NULL) {
715:           comm = this->comm();
716:           rank = this->commRank();
717:         } else {
718:           MPI_Comm_rank(comm, &rank);
719:         }
720:         if (name == "") {
721:           if(rank == 0) {
722:             txt << "viewing a ConeSection" << std::endl;
723:           }
724:         } else {
725:           if(rank == 0) {
726:             txt << "viewing ConeSection '" << name << "'" << std::endl;
727:           }
728:         }
729:         PetscSynchronizedPrintf(comm, txt.str().c_str());
730:         PetscSynchronizedFlush(comm);
731:       };
732:     };

734:     template<typename Bundle_, typename Sieve_>
735:     class SupportSizeSection : public ALE::ParallelObject {
736:     public:
737:       typedef Bundle_                          bundle_type;
738:       typedef Sieve_                           sieve_type;
739:       typedef typename sieve_type::source_type point_type;
740:       typedef typename sieve_type::target_type value_type;
741:     protected:
742:       Obj<bundle_type> _bundle;
743:       Obj<sieve_type>  _sieve;
744:       value_type       _size;
745:       int              _minDepth;
746:     public:
747:       SupportSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumDepth = 0) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _bundle(bundle), _sieve(sieve), _minDepth(minimumDepth) {};
748:       virtual ~SupportSizeSection() {};
749:     public:
750:       bool hasPoint(const point_type& point) {return true;};
751:       const value_type *restrictPoint(const point_type& p) {
752:         if ((this->_minDepth == 0) || (this->_bundle->depth(p) >= this->_minDepth)) {
753:           this->_size = this->_sieve->support(p)->size();
754:         } else {
755:           this->_size = 0;
756:         }
757:         return &this->_size;
758:       };
759:     public:
760:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
761:         ostringstream txt;
762:         int rank;

764:         if (comm == MPI_COMM_NULL) {
765:           comm = this->comm();
766:           rank = this->commRank();
767:         } else {
768:           MPI_Comm_rank(comm, &rank);
769:         }
770:         if (name == "") {
771:           if(rank == 0) {
772:             txt << "viewing a SupportSizeSection" << std::endl;
773:           }
774:         } else {
775:           if(rank == 0) {
776:             txt << "viewing SupportSizeSection '" << name << "'" << std::endl;
777:           }
778:         }
779:         PetscSynchronizedPrintf(comm, txt.str().c_str());
780:         PetscSynchronizedFlush(comm);
781:       };
782:     };

784:     template<typename Sieve_>
785:     class SupportSection : public ALE::ParallelObject {
786:     public:
787:       typedef Sieve_                           sieve_type;
788:       typedef typename sieve_type::source_type point_type;
789:       typedef typename sieve_type::target_type value_type;
790:       typedef PartitionDomain<sieve_type>      chart_type;
791:     protected:
792:       Obj<sieve_type> _sieve;
793:       int             _supportSize;
794:       value_type     *_support;
795:       chart_type      _domain;
796:       void ensureSupport(const int size) {
797:         if (size > this->_supportSize) {
798:           if (this->_support) delete [] this->_support;
799:           this->_supportSize = size;
800:           this->_support     = new value_type[this->_supportSize];
801:         }
802:       };
803:     public:
804:       SupportSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _supportSize(-1), _support(NULL) {};
805:       virtual ~SupportSection() {if (this->_support) delete [] this->_support;};
806:     public:
807:       const chart_type& getChart() {return this->_domain;};
808:       bool hasPoint(const point_type& point) {return true;};
809:       const value_type *restrictPoint(const point_type& p) {
810:         const Obj<typename sieve_type::traits::supportSequence>& support = this->_sieve->support(p);
811:         int s = 0;

813:         this->ensureSupport(support->size());
814:         for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
815:           this->_support[s++] = *s_iter;
816:         }
817:         return this->_support;
818:       };
819:     public:
820:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
821:         ostringstream txt;
822:         int rank;

824:         if (comm == MPI_COMM_NULL) {
825:           comm = this->comm();
826:           rank = this->commRank();
827:         } else {
828:           MPI_Comm_rank(comm, &rank);
829:         }
830:         if (name == "") {
831:           if(rank == 0) {
832:             txt << "viewing a SupportSection" << std::endl;
833:           }
834:         } else {
835:           if(rank == 0) {
836:             txt << "viewing SupportSection '" << name << "'" << std::endl;
837:           }
838:         }
839:         PetscSynchronizedPrintf(comm, txt.str().c_str());
840:         PetscSynchronizedFlush(comm);
841:       };
842:     };

844:     template<typename Sieve_, typename Section_>
845:     class ArrowSection : public ALE::ParallelObject {
846:     public:
847:       typedef Sieve_                            sieve_type;
848:       typedef Section_                          section_type;
849:       typedef typename sieve_type::target_type  point_type;
850:       typedef typename section_type::point_type arrow_type;
851:       typedef typename section_type::value_type value_type;
852:     protected:
853:       Obj<sieve_type>   _sieve;
854:       Obj<section_type> _section;
855:       int               _coneSize;
856:       value_type       *_cone;
857:       void ensureCone(const int size) {
858:         if (size > this->_coneSize) {
859:           if (this->_cone) delete [] this->_cone;
860:           this->_coneSize = size;
861:           this->_cone     = new value_type[this->_coneSize];
862:         }
863:       };
864:     public:
865:       ArrowSection(const Obj<sieve_type>& sieve, const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _section(section), _coneSize(-1), _cone(NULL) {};
866:       virtual ~ArrowSection() {if (this->_cone) delete [] this->_cone;};
867:     public:
868:       bool hasPoint(const point_type& point) {return this->_sieve->baseContains(point);};
869:       const value_type *restrictPoint(const point_type& p) {
870:         const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
871:         int c = -1;

873:         this->ensureCone(cone->size());
874:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
875:           this->_cone[++c] = this->_section->restrictPoint(arrow_type(*c_iter, p))[0];
876:         }
877:         return this->_cone;
878:       };
879:     public:
880:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
881:         ostringstream txt;
882:         int rank;

884:         if (comm == MPI_COMM_NULL) {
885:           comm = this->comm();
886:           rank = this->commRank();
887:         } else {
888:           MPI_Comm_rank(comm, &rank);
889:         }
890:         if (name == "") {
891:           if(rank == 0) {
892:             txt << "viewing a ConeSection" << std::endl;
893:           }
894:         } else {
895:           if(rank == 0) {
896:             txt << "viewing ConeSection '" << name << "'" << std::endl;
897:           }
898:         }
899:         PetscSynchronizedPrintf(comm, txt.str().c_str());
900:         PetscSynchronizedFlush(comm);
901:       };
902:     };
903:   }
904: }
905: #endif