Actual source code: ISieve.hh

  1: #ifndef included_ALE_ISieve_hh
  2: #define included_ALE_ISieve_hh

  4: #ifndef  included_ALE_hh
  5: #include <ALE.hh>
  6: #endif

  8: #include <fstream>

 10: //#define IMESH_NEW_LABELS

 12: namespace ALE {
 13:   template<typename Point>
 14:   class OrientedPoint : public std::pair<Point, int> {
 15:   public:
 16:     OrientedPoint(const int o) : std::pair<Point, int>(o, o) {};
 17:     ~OrientedPoint() {};
 18:     friend std::ostream& operator<<(std::ostream& stream, const OrientedPoint& point) {
 19:       stream << "(" << point.first << ", " << point.second << ")";
 20:       return stream;
 21:     };
 22:   };

 24:   template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
 25:   class Interval {
 26:   public:
 27:     typedef Point_            point_type;
 28:     typedef Alloc_            alloc_type;
 29:   public:
 30:     class const_iterator {
 31:     protected:
 32:       point_type _p;
 33:     public:
 34:       const_iterator(const point_type p): _p(p) {};
 35:       ~const_iterator() {};
 36:     public:
 37:       const_iterator& operator=(const const_iterator& iter) {this->_p = iter._p; return *this;};
 38:       bool operator==(const const_iterator& iter) const {return this->_p == iter._p;};
 39:       bool operator!=(const const_iterator& iter) const {return this->_p != iter._p;};
 40:       const_iterator& operator++() {++this->_p; return *this;}
 41:       const_iterator& operator++(int) {
 42:         const_iterator tmp(*this);
 43:         ++(*this);
 44:         return tmp;
 45:       };
 46:       const_iterator& operator--() {--this->_p; return *this;}
 47:       const_iterator& operator--(int) {
 48:         const_iterator tmp(*this);
 49:         --(*this);
 50:         return tmp;
 51:       };
 52:       point_type operator*() const {return this->_p;};
 53:     };
 54:   protected:
 55:     point_type _min, _max;
 56:   public:
 57:     Interval(): _min(point_type()), _max(point_type()) {};
 58:     Interval(const point_type& min, const point_type& max): _min(min), _max(max) {};
 59:     Interval(const Interval& interval): _min(interval.min()), _max(interval.max()) {};
 60:     template<typename Iterator>
 61:     Interval(Iterator& iterator) {
 62:       this->_min = *std::min_element(iterator.begin(), iterator.end());
 63:       this->_max = (*std::max_element(iterator.begin(), iterator.end()))+1;
 64:     }
 65:   public:
 66:     Interval& operator=(const Interval& interval) {_min = interval.min(); _max = interval.max(); return *this;}
 67:     friend std::ostream& operator<<(std::ostream& stream, const Interval& interval) {
 68:       stream << "(" << interval.min() << ", " << interval.max() << ")";
 69:       return stream;
 70:     }
 71:   public:
 72:     const_iterator begin() const {return const_iterator(this->_min);};
 73:     const_iterator end() const {return const_iterator(this->_max);};
 74:     size_t size() const {return this->_max - this->_min;};
 75:     size_t count(const point_type& p) const {return ((p >= _min) && (p < _max)) ? 1 : 0;};
 76:     point_type min() const {return this->_min;};
 77:     point_type max() const {return this->_max;};
 78:     bool hasPoint(const point_type& point) const {
 79:       if (point < this->_min || point >= this->_max) return false;
 80:       return true;
 81:     };
 82:     void checkPoint(const point_type& point) const {
 83:       if (point < this->_min || point >= this->_max) {
 84:         ostringstream msg;
 85:         msg << "Invalid point " << point << " not in " << *this << std::endl;
 86:         throw ALE::Exception(msg.str().c_str());
 87:       }
 88:     };
 89:   };

 91:   template<typename Source_, typename Target_>
 92:   struct SimpleArrow {
 93:     typedef Source_ source_type;
 94:     typedef Target_ target_type;
 95:     const source_type source;
 96:     const target_type target;
 97:     SimpleArrow(const source_type& s, const target_type& t) : source(s), target(t) {};
 98:     template<typename OtherSource_, typename OtherTarget_>
 99:     struct rebind {
100:       typedef SimpleArrow<OtherSource_, OtherTarget_> other;
101:     };
102:     struct flip {
103:       typedef SimpleArrow<target_type, source_type> other;
104:       other arrow(const SimpleArrow& a) {return type(a.target, a.source);};
105:     };
106:     friend bool operator<(const SimpleArrow& x, const SimpleArrow& y) {
107:       return ((x.source < y.source) || ((x.source == y.source) && (x.target < y.target)));
108:     };
109:     friend std::ostream& operator<<(std::ostream& os, const SimpleArrow& a) {
110:       os << a.source << " ----> " << a.target;
111:       return os;
112:     }
113:   };

115:   namespace ISieveVisitor {
116:     template<typename Sieve>
117:     class NullVisitor {
118:     public:
119:       inline void visitArrow(const typename Sieve::arrow_type&) {};
120:       inline void visitPoint(const typename Sieve::point_type&) {};
121:       inline void visitArrow(const typename Sieve::arrow_type&, const int orientation) {};
122:       inline void visitPoint(const typename Sieve::point_type&, const int orientation) {};
123:     };
124:     class PrintVisitor {
125:     protected:
126:       ostringstream& os;
127:       const int      rank;
128:     public:
129:       PrintVisitor(ostringstream& s, const int rank = 0) : os(s), rank(rank) {};
130:       template<typename Arrow>
131:       inline void visitArrow(const Arrow& arrow) const {
132:         this->os << "["<<this->rank<<"]: " << arrow << std::endl;
133:       }
134:       template<typename Point>
135:       inline void visitPoint(const Point&) const {}
136:     };
137:     class ReversePrintVisitor : public PrintVisitor {
138:     public:
139:       ReversePrintVisitor(ostringstream& s, const int rank) : PrintVisitor(s, rank) {};
140:       template<typename Arrow>
141:       inline void visitArrow(const Arrow& arrow) const {
142:         this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << std::endl;
143:       }
144:       template<typename Arrow>
145:       inline void visitArrow(const Arrow& arrow, const int orientation) const {
146:         this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << ": " << orientation << std::endl;
147:       }
148:       template<typename Point>
149:       inline void visitPoint(const Point&) const {}
150:       template<typename Point>
151:       inline void visitPoint(const Point&, const int) const {}
152:     };
153:     template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
154:     class PointRetriever {
155:     public:
156:       typedef typename Sieve::point_type point_type;
157:       typedef typename Sieve::arrow_type arrow_type;
158:       typedef std::pair<point_type,int>  oriented_point_type;
159:     protected:
160:       const bool           unique;
161:       size_t               i, o;
162:       size_t               skip, limit;
163:       Visitor             *visitor;
164:       size_t               size;
165:       point_type          *points;
166:       oriented_point_type *oPoints;
167:     protected:
168:       inline virtual bool accept(const point_type& point) {return true;};
169:     public:
170:       PointRetriever() : unique(false), i(0), o(0), skip(0), limit(0) {
171:         this->size    = 0;
172:         this->points  = NULL;
173:         this->oPoints = NULL;
174:       };
175:       PointRetriever(const size_t size, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0) {
176:         static Visitor nV;
177:         this->visitor = &nV;
178:         this->points  = NULL;
179:         this->oPoints = NULL;
180:         this->setSize(size);
181:       };
182:       PointRetriever(const size_t size, Visitor& v, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0), visitor(&v) {
183:         this->points  = NULL;
184:         this->oPoints = NULL;
185:         this->setSize(size);
186:       };
187:       virtual ~PointRetriever() {
188:         delete [] this->points;
189:         delete [] this->oPoints;
190:         this->points  = NULL;
191:         this->oPoints = NULL;
192:       };
193:       inline void visitArrow(const arrow_type& arrow) {
194:         this->visitor->visitArrow(arrow);
195:       };
196:       inline void visitArrow(const arrow_type& arrow, const int orientation) {
197:         this->visitor->visitArrow(arrow, orientation);
198:       };
199:       inline void visitPoint(const point_type& point) {
200:         if (i >= size) {
201:           ostringstream msg;
202:           msg << "Too many points (>" << size << ")for PointRetriever visitor";
203:           throw ALE::Exception(msg.str().c_str());
204:         }
205:         if (this->accept(point)) {
206:           if (this->unique) {
207:             size_t p;
208:             for(p = 0; p < i; ++p) {if (points[p] == point) break;}
209:             if (p != i) return;
210:           }
211:           if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
212:           points[i++] = point;
213:           this->visitor->visitPoint(point);
214:         }
215:       };
216:       inline void visitPoint(const point_type& point, const int orientation) {
217:         if (o >= size) {
218:           ostringstream msg;
219:           msg << "Too many ordered points (>" << size << ")for PointRetriever visitor";
220:           throw ALE::Exception(msg.str().c_str());
221:         }
222:         if (this->accept(point)) {
223:           if (this->unique) {
224:             size_t p;
225:             for(p = 0; p < i; ++p) {if (points[p] == point) break;}
226:             if (p != i) return;
227:           }
228:           if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
229:           points[i++]  = point;
230:           oPoints[o++] = oriented_point_type(point, orientation);
231:           this->visitor->visitPoint(point, orientation);
232:         }
233:       };
234:     public:
235:       size_t                     getSize() const {return this->i;}
236:       const point_type          *getPoints() const {return this->points;}
237:       size_t                     getOrientedSize() const {return this->o;}
238:       const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
239:       void clear() {this->i = this->o = 0;}
240:       void setSize(const size_t s) {
241:         if (this->points) {
242:           delete [] this->points;
243:           delete [] this->oPoints;
244:         }
245:         this->size    = s;
246:         this->points  = new point_type[this->size];
247:         this->oPoints = new oriented_point_type[this->size];
248:       }
249:       void setSkip(size_t s) {this->skip = s;};
250:       void setLimit(size_t l) {this->limit = l;};
251:     };
252:     template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
253:     class NConeRetriever : public PointRetriever<Sieve,Visitor> {
254:     public:
255:       typedef PointRetriever<Sieve,Visitor>           base_type;
256:       typedef typename Sieve::point_type              point_type;
257:       typedef typename Sieve::arrow_type              arrow_type;
258:       typedef typename base_type::oriented_point_type oriented_point_type;
259:     protected:
260:       const Sieve& sieve;
261:     protected:
262:       inline virtual bool accept(const point_type& point) {
263:         if (!this->sieve.getConeSize(point))
264:           return true;
265:         return false;
266:       };
267:     public:
268:       NConeRetriever(const Sieve& s, const size_t size) : PointRetriever<Sieve,Visitor>(size, true), sieve(s) {};
269:       NConeRetriever(const Sieve& s, const size_t size, Visitor& v) : PointRetriever<Sieve,Visitor>(size, v, true), sieve(s) {};
270:       virtual ~NConeRetriever() {};
271:     };
272:     template<typename Mesh, typename Visitor = NullVisitor<typename Mesh::sieve_type> >
273:     class MeshNConeRetriever : public PointRetriever<typename Mesh::sieve_type,Visitor> {
274:     public:
275:       typedef typename Mesh::Sieve                    Sieve;
276:       typedef PointRetriever<Sieve,Visitor>           base_type;
277:       typedef typename Sieve::point_type              point_type;
278:       typedef typename Sieve::arrow_type              arrow_type;
279:       typedef typename base_type::oriented_point_type oriented_point_type;
280:     protected:
281:       const Mesh& mesh;
282:       const int   depth;
283:     protected:
284:       inline virtual bool accept(const point_type& point) {
285:         if (this->mesh.depth(point) == this->depth)
286:           return true;
287:         return false;
288:       };
289:     public:
290:       MeshNConeRetriever(const Mesh& m, const int depth, const size_t size) : PointRetriever<typename Mesh::Sieve,Visitor>(size, true), mesh(m), depth(depth) {};
291:       MeshNConeRetriever(const Mesh& m, const int depth, const size_t size, Visitor& v) : PointRetriever<typename Mesh::Sieve,Visitor>(size, v, true), mesh(m), depth(depth) {};
292:       virtual ~MeshNConeRetriever() {};
293:     };
294:     template<typename Sieve, typename Set, typename Renumbering>
295:     class FilteredPointRetriever {
296:     public:
297:       typedef typename Sieve::point_type point_type;
298:       typedef typename Sieve::arrow_type arrow_type;
299:       typedef std::pair<point_type,int>  oriented_point_type;
300:     protected:
301:       const Set&   pointSet;
302:       Renumbering& renumbering;
303:       const size_t size;
304:       size_t       i, o;
305:       bool         renumber;
306:       point_type  *points;
307:       oriented_point_type *oPoints;
308:     public:
309:       FilteredPointRetriever(const Set& s, Renumbering& r, const size_t size) : pointSet(s), renumbering(r), size(size), i(0), o(0), renumber(true) {
310:         this->points  = new point_type[this->size];
311:         this->oPoints = new oriented_point_type[this->size];
312:       };
313:       ~FilteredPointRetriever() {
314:         delete [] this->points;
315:         delete [] this->oPoints;
316:       };
317:       inline void visitArrow(const arrow_type& arrow) {};
318:       inline void visitPoint(const point_type& point) {
319:         if (i >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
320:         if (this->pointSet.find(point) == this->pointSet.end()) return;
321:         if (renumber) {
322:           points[i++] = this->renumbering[point];
323:         } else {
324:           points[i++] = point;
325:         }
326:       };
327:       inline void visitArrow(const arrow_type& arrow, const int orientation) {};
328:       inline void visitPoint(const point_type& point, const int orientation) {
329:         if (o >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
330:         if (this->pointSet.find(point) == this->pointSet.end()) return;
331:         if (renumber) {
332:           points[i++]  = this->renumbering[point];
333:           oPoints[o++] = oriented_point_type(this->renumbering[point], orientation);
334:         } else {
335:           points[i++]  = point;
336:           oPoints[o++] = oriented_point_type(point, orientation);
337:         }
338:       };
339:     public:
340:       size_t            getSize() const {return this->i;}
341:       const point_type *getPoints() const {return this->points;}
342:       size_t            getOrientedSize() const {return this->o;}
343:       const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
344:       void clear() {this->i = 0; this->o = 0;}
345:       void useRenumbering(const bool renumber) {this->renumber = renumber;}
346:     };
347:     template<typename Sieve, int size, typename Visitor = NullVisitor<Sieve> >
348:     class ArrowRetriever {
349:     public:
350:       typedef typename Sieve::point_type point_type;
351:       typedef typename Sieve::arrow_type arrow_type;
352:       typedef std::pair<arrow_type,int>  oriented_arrow_type;
353:     protected:
354:       arrow_type          arrows[size];
355:       oriented_arrow_type oArrows[size];
356:       size_t              i, o;
357:       Visitor            *visitor;
358:     public:
359:       ArrowRetriever() : i(0), o(0) {
360:         static Visitor nV;
361:         this->visitor = &nV;
362:       };
363:       ArrowRetriever(Visitor& v) : i(0), o(0), visitor(&v) {};
364:       inline void visitArrow(const typename Sieve::arrow_type& arrow) {
365:         if (i >= size) throw ALE::Exception("Too many arrows for visitor");
366:         arrows[i++] = arrow;
367:         this->visitor->visitArrow(arrow);
368:       };
369:       inline void visitArrow(const typename Sieve::arrow_type& arrow, const int orientation) {
370:         if (o >= size) throw ALE::Exception("Too many arrows for visitor");
371:         oArrows[o++] = oriented_arrow_type(arrow, orientation);
372:         this->visitor->visitArrow(arrow, orientation);
373:       };
374:       inline void visitPoint(const point_type& point) {
375:         this->visitor->visitPoint(point);
376:       };
377:       inline void visitPoint(const point_type& point, const int orientation) {
378:         this->visitor->visitPoint(point, orientation);
379:       };
380:     public:
381:       size_t                     getSize() const {return this->i;}
382:       const point_type          *getArrows() const {return this->arrows;}
383:       size_t                     getOrientedSize() const {return this->o;}
384:       const oriented_arrow_type *getOrientedArrows() const {return this->oArrows;}
385:       void clear() {this->i = this->o = 0;}
386:     };
387:     template<typename Sieve, typename Visitor>
388:     class ConeVisitor {
389:     protected:
390:       const Sieve& sieve;
391:       Visitor&     visitor;
392:       bool         useSource;
393:     public:
394:       ConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
395:       inline void visitPoint(const typename Sieve::point_type& point) {
396:         this->sieve.cone(point, visitor);
397:       };
398:       inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
399:     };
400:     template<typename Sieve, typename Visitor>
401:     class OrientedConeVisitor {
402:     protected:
403:       const Sieve& sieve;
404:       Visitor&     visitor;
405:       bool         useSource;
406:     public:
407:       OrientedConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
408:       inline void visitPoint(const typename Sieve::point_type& point) {
409:         this->sieve.orientedCone(point, visitor);
410:       };
411:       inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
412:     };
413:     template<typename Sieve, typename Visitor>
414:     class SupportVisitor {
415:     protected:
416:       const Sieve& sieve;
417:       Visitor&     visitor;
418:       bool         useSource;
419:     public:
420:       SupportVisitor(const Sieve& s, Visitor& v, bool useSource = true) : sieve(s), visitor(v), useSource(useSource) {};
421:       inline void visitPoint(const typename Sieve::point_type& point) {
422:         this->sieve.support(point, visitor);
423:       };
424:       inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
425:     };
426:     template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
427:     class TransitiveClosureVisitor {
428:     public:
429:       typedef Visitor visitor_type;
430:     protected:
431:       const Sieve& sieve;
432:       Visitor&     visitor;
433:       bool         isCone;
434:       std::set<typename Sieve::point_type> seen;
435:     public:
436:       TransitiveClosureVisitor(const Sieve& s, Visitor& v) : sieve(s), visitor(v), isCone(true) {};
437:       inline void visitPoint(const typename Sieve::point_type& point) const {};
438:       inline void visitArrow(const typename Sieve::arrow_type& arrow) {
439:         if (this->isCone) {
440:           if (this->seen.find(arrow.target) == this->seen.end()) {
441:             this->seen.insert(arrow.target);
442:             this->visitor.visitPoint(arrow.target);
443:           }
444:           this->visitor.visitArrow(arrow);
445:           if (this->seen.find(arrow.source) == this->seen.end()) {
446:             if (this->sieve.getConeSize(arrow.source)) {
447:               this->sieve.cone(arrow.source, *this);
448:             } else {
449:               this->seen.insert(arrow.source);
450:               this->visitor.visitPoint(arrow.source);
451:             }
452:           }
453:         } else {
454:           if (this->seen.find(arrow.source) == this->seen.end()) {
455:             this->seen.insert(arrow.source);
456:             this->visitor.visitPoint(arrow.source);
457:           }
458:           this->visitor.visitArrow(arrow);
459:           if (this->seen.find(arrow.target) == this->seen.end()) {
460:             if (this->sieve.getSupportSize(arrow.target)) {
461:               this->sieve.support(arrow.target, *this);
462:             } else {
463:               this->seen.insert(arrow.target);
464:               this->visitor.visitPoint(arrow.target);
465:             }
466:           }
467:         }
468:       };
469:     public:
470:       bool getIsCone() const {return this->isCone;};
471:       void setIsCone(const bool isCone) {this->isCone = isCone;};
472:       const std::set<typename Sieve::point_type>& getPoints() const {return this->seen;};
473:       void clear() {this->seen.clear();};
474:     };
475:     template<typename Sieve, typename Section>
476:     class SizeVisitor {
477:     protected:
478:       const Section& section;
479:       int            size;
480:     public:
481:       SizeVisitor(const Section& s) : section(s), size(0) {};
482:       inline void visitPoint(const typename Sieve::point_type& point) {
483:         this->size += section.getConstrainedFiberDimension(point);
484:       };
485:       inline void visitArrow(const typename Sieve::arrow_type&) {};
486:     public:
487:       int getSize() {return this->size;};
488:     };
489:     template<typename Sieve, typename Section>
490:     class SizeWithBCVisitor {
491:     protected:
492:       const Section& section;
493:       int            size;
494:     public:
495:       SizeWithBCVisitor(const Section& s) : section(s), size(0) {};
496:       inline void visitPoint(const typename Sieve::point_type& point) {
497:         this->size += section.getFiberDimension(point);
498:       };
499:       inline void visitArrow(const typename Sieve::arrow_type&) {};
500:     public:
501:       int getSize() {return this->size;};
502:     };
503:     template<typename Section>
504:     class RestrictVisitor {
505:     public:
506:       typedef typename Section::value_type value_type;
507:     protected:
508:       const Section& section;
509:       int            size;
510:       int            i;
511:       value_type    *values;
512:       bool           allocated;
513:     public:
514:       RestrictVisitor(const Section& s, const int size) : section(s), size(size), i(0) {
515:         this->values    = new value_type[this->size];
516:         this->allocated = true;
517:       };
518:       RestrictVisitor(const Section& s, const int size, value_type *values) : section(s), size(size), i(0) {
519:         this->values    = values;
520:         this->allocated = false;
521:       };
522:       ~RestrictVisitor() {if (this->allocated) {delete [] this->values;}};
523:       template<typename Point>
524:       inline void visitPoint(const Point& point, const int orientation) {
525:         const int         dim = section.getFiberDimension(point);
526:         if (i+dim > size) {throw ALE::Exception("Too many values for RestrictVisitor.");}
527:         const value_type *v   = section.restrictPoint(point);

529:         if (orientation >= 0) {
530:           for(int d = 0; d < dim; ++d, ++i) {
531:             this->values[i] = v[d];
532:           }
533:         } else {
534:           for(int d = dim-1; d >= 0; --d, ++i) {
535:             this->values[i] = v[d];
536:           }
537:         }
538:       }
539:       template<typename Arrow>
540:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
541:     public:
542:       const value_type *getValues() const {return this->values;};
543:       int  getSize() const {return this->i;};
544:       int  getMaxSize() const {return this->size;};
545:       void ensureSize(const int size) {
546:         this->clear();
547:         if (size > this->size) {
548:           this->size = size;
549:           if (this->allocated) {delete [] this->values;}
550:           this->values = new value_type[this->size];
551:           this->allocated = true;
552:         }
553:       };
554:       void clear() {this->i = 0;};
555:     };
556:     template<typename Section>
557:     class UpdateVisitor {
558:     public:
559:       typedef typename Section::value_type value_type;
560:     protected:
561:       Section&          section;
562:       const value_type *values;
563:       int               i;
564:     public:
565:       UpdateVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
566:       template<typename Point>
567:       inline void visitPoint(const Point& point, const int orientation) {
568:         const int dim = section.getFiberDimension(point);
569:         this->section.updatePoint(point, &this->values[this->i], orientation);
570:         this->i += dim;
571:       }
572:       template<typename Arrow>
573:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
574:       void clear() {this->i = 0;};
575:     };
576:     template<typename Section>
577:     class UpdateAllVisitor {
578:     public:
579:       typedef typename Section::value_type value_type;
580:     protected:
581:       Section&          section;
582:       const value_type *values;
583:       int               i;
584:     public:
585:       UpdateAllVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
586:       template<typename Point>
587:       inline void visitPoint(const Point& point, const int orientation) {
588:         const int dim = section.getFiberDimension(point);
589:         this->section.updatePointAll(point, &this->values[this->i], orientation);
590:         this->i += dim;
591:       }
592:       template<typename Arrow>
593:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
594:       void clear() {this->i = 0;};
595:     };
596:     template<typename Section>
597:     class UpdateAddVisitor {
598:     public:
599:       typedef typename Section::value_type value_type;
600:     protected:
601:       Section&          section;
602:       const value_type *values;
603:       int               i;
604:     public:
605:       UpdateAddVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
606:       template<typename Point>
607:       inline void visitPoint(const Point& point, const int orientation) {
608:         const int dim = section.getFiberDimension(point);
609:         this->section.updateAddPoint(point, &this->values[this->i], orientation);
610:         this->i += dim;
611:       }
612:       template<typename Arrow>
613:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
614:       void clear() {this->i = 0;};
615:     };
616:     template<typename Section, typename Order, typename Value>
617:     class IndicesVisitor {
618:     public:
619:       typedef Value                        value_type;
620:       typedef typename Section::point_type point_type;
621:     protected:
622:       const Section& section;
623:       // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
624:       Order&         order;
625:       int            size;
626:       int            i, p;
627:       // If false, constrained indices are returned as negative values. Otherwise, they are omitted
628:       bool           freeOnly;
629:       // If true, it allows space for constrained variables (even if the indices are not returned) Wierd
630:       bool           skipConstraints;
631:       value_type    *values;
632:       bool           allocated;
633:       point_type    *points;
634:       bool           allocatedPoints;
635:     protected:
636:       void getUnconstrainedIndices(const point_type& p, const int orientation) {
637:         if (i+section.getFiberDimension(p) > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
638:         if (orientation >= 0) {
639:           const int start = this->order.getIndex(p);
640:           const int end   = start + section.getFiberDimension(p);

642:           for(int j = start; j < end; ++j, ++i) {
643:             this->values[i] = j;
644:           }
645:         } else if (!section.getNumSpaces()) {
646:           const int start = this->order.getIndex(p);
647:           const int end   = start + section.getFiberDimension(p);

649:           for(int j = end-1; j >= start; --j, ++i) {
650:             this->values[i] = j;
651:           }
652:         } else {
653:           const int numSpaces = section.getNumSpaces();
654:           int       start     = this->order.getIndex(p);

656:           for(int space = 0; space < numSpaces; ++space) {
657:             const int end = start + section.getFiberDimension(p, space);

659:             for(int j = end-1; j >= start; --j, ++i) {
660:               this->values[i] = j;
661:             }
662:             start = end;
663:           }
664:         }
665:       };
666:       void getConstrainedIndices(const point_type& p, const int orientation) {
667:         const int cDim = this->section.getConstraintDimension(p);
668:         if (i+cDim > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
669:         typedef typename Section::bc_type::value_type index_type;
670:         const index_type *cDof  = this->section.getConstraintDof(p);
671:         const int         start = this->order.getIndex(p);

673:         if (orientation >= 0) {
674:           const int dim = this->section.getFiberDimension(p);

676:           for(int j = start, cInd = 0, k = 0; k < dim; ++k) {
677:             if ((cInd < cDim) && (k == cDof[cInd])) {
678:               if (!freeOnly) values[i++] = -(k+1);
679:               if (skipConstraints) ++j;
680:               ++cInd;
681:             } else {
682:               values[i++] = j++;
683:             }
684:           }
685:         } else {
686:           int offset  = 0;
687:           int cOffset = 0;
688:           int k       = -1;

690:           for(int space = 0; space < section.getNumSpaces(); ++space) {
691:             const int  dim = this->section.getFiberDimension(p, space);
692:             const int tDim = this->section.getConstrainedFiberDimension(p, space);
693:             int       cInd = (dim - tDim)-1;

695:             k += dim;
696:             for(int l = 0, j = start+tDim+offset; l < dim; ++l, --k) {
697:               if ((cInd >= 0) && (k == cDof[cInd+cOffset])) {
698:                 if (!freeOnly) values[i++] = -(offset+l+1);
699:                 if (skipConstraints) --j;
700:                 --cInd;
701:               } else {
702:                 values[i++] = --j;
703:               }
704:             }
705:             k       += dim;
706:             offset  += dim;
707:             cOffset += dim - tDim;
708:           }
709:         }
710:       };
711:     public:
712:       IndicesVisitor(const Section& s, Order& o, const int size, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
713:         this->values    = new value_type[this->size];
714:         this->allocated = true;
715:         if (unique) {
716:           this->points          = new point_type[this->size];
717:           this->allocatedPoints = true;
718:         } else {
719:           this->points          = NULL;
720:           this->allocatedPoints = false;
721:         }
722:       };
723:       IndicesVisitor(const Section& s, Order& o, const int size, value_type *values, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
724:         this->values    = values;
725:         this->allocated = false;
726:         if (unique) {
727:           this->points          = new point_type[this->size];
728:           this->allocatedPoints = true;
729:         } else {
730:           this->points          = NULL;
731:           this->allocatedPoints = false;
732:         }
733:       };
734:       ~IndicesVisitor() {
735:         if (this->allocated) {delete [] this->values;}
736:         if (this->allocatedPoints) {delete [] this->points;}
737:       };
738:       inline void visitPoint(const point_type& point, const int orientation) {
739:         if (p >= size) {
740:           ostringstream msg;
741:           msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
742:           throw ALE::Exception(msg.str().c_str());
743:         }
744:         if (points) {
745:           int pp;
746:           for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
747:           if (pp != p) return;
748:           points[p++] = point;
749:         }
750:         const int cDim = this->section.getConstraintDimension(point);

752:         if (!cDim) {
753:           this->getUnconstrainedIndices(point, orientation);
754:         } else {
755:           this->getConstrainedIndices(point, orientation);
756:         }
757:       }
758:       template<typename Arrow>
759:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
760:     public:
761:       const value_type *getValues() const {return this->values;};
762:       int  getSize() const {return this->i;};
763:       int  getMaxSize() const {return this->size;};
764:       void ensureSize(const int size) {
765:         this->clear();
766:         if (size > this->size) {
767:           this->size = size;
768:           if (this->allocated) {delete [] this->values;}
769:           this->values = new value_type[this->size];
770:           this->allocated = true;
771:           if (this->allocatedPoints) {delete [] this->points;}
772:           this->points = new point_type[this->size];
773:           this->allocatedPoints = true;
774:         }
775:       };
776:       void clear() {this->i = 0; this->p = 0;};
777:     };
778:     template<typename Sieve, typename Label>
779:     class MarkVisitor {
780:     protected:
781:       Label& label;
782:       int    marker;
783:     public:
784:       MarkVisitor(Label& l, const int marker) : label(l), marker(marker) {};
785:       inline void visitPoint(const typename Sieve::point_type& point) {
786:         this->label.setCone(this->marker, point);
787:       };
788:       inline void visitArrow(const typename Sieve::arrow_type&) {};
789:     };
790:   };

792:   template<typename Sieve>
793:   class ISieveTraversal {
794:   public:
795:     typedef typename Sieve::point_type point_type;
796:   public:
797:     template<typename Visitor>
798:     static void orientedClosure(const Sieve& sieve, const point_type& p, Visitor& v) {
799:       typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
800:       Retriever cV[2] = {Retriever(200,v), Retriever(200,v)};
801:       int       c     = 0;

803:       v.visitPoint(p, 0);
804:       // Cone is guarateed to be ordered correctly
805:       sieve.orientedCone(p, cV[c]);

807:       while(cV[c].getOrientedSize()) {
808:         const typename Retriever::oriented_point_type *cone     = cV[c].getOrientedPoints();
809:         const int                                      coneSize = cV[c].getOrientedSize();
810:         c = 1 - c;

812:         for(int p = 0; p < coneSize; ++p) {
813:           const typename Retriever::point_type& point     = cone[p].first;
814:           int                                   pO        = cone[p].second == 0 ? 1 : cone[p].second;
815:           const int                             pConeSize = sieve.getConeSize(point);

817:           if (pO < 0) {
818:             if (pO == -pConeSize) {
819:               sieve.orientedReverseCone(point, cV[c]);
820:             } else {
821:               const int numSkip = sieve.getConeSize(point) + pO;

823:               cV[c].setSkip(cV[c].getSize()+numSkip);
824:               cV[c].setLimit(cV[c].getSize()+pConeSize);
825:               sieve.orientedReverseCone(point, cV[c]);
826:               sieve.orientedReverseCone(point, cV[c]);
827:               cV[c].setSkip(0);
828:               cV[c].setLimit(0);
829:             }
830:           } else {
831:             if (pO == 1) {
832:               sieve.orientedCone(point, cV[c]);
833:             } else {
834:               const int numSkip = pO-1;

836:               cV[c].setSkip(cV[c].getSize()+numSkip);
837:               cV[c].setLimit(cV[c].getSize()+pConeSize);
838:               sieve.orientedCone(point, cV[c]);
839:               sieve.orientedCone(point, cV[c]);
840:               cV[c].setSkip(0);
841:               cV[c].setLimit(0);
842:             }
843:           }
844:         }
845:         cV[1-c].clear();
846:       }
847: #if 0
848:       // These contain arrows paired with orientations from the \emph{previous} arrow
849:       Obj<orientedArrowArray> cone    = new orientedArrowArray();
850:       Obj<orientedArrowArray> base    = new orientedArrowArray();
851:       coneSet                 seen;

853:       for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
854:         const arrow_type arrow(*c_iter, p);

856:         cone->push_back(oriented_arrow_type(arrow, 1)); // Notice the orientation of leaf cells is always 1
857:         closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0])); // However, we return the actual arrow orientation
858:       }
859:       for(int i = 1; i < depth; ++i) {
860:         Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;

862:         cone->clear();
863:         for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
864:           const arrow_type&                                     arrow = b_iter->first; // This is an arrow considered in the previous round
865:           const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source); // We are going to get the cone of this guy
866:           typename arrow_section_type::value_type               o     = orientation->restrictPoint(arrow)[0]; // The orientation of arrow, which is our pO

868:           if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
869:             o = -(o+1);
870:           }
871:           if (o < 0) {
872:             const int size = pCone->size();

874:             if (o == -size) {
875:               for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
876:                 if (seen.find(*c_iter) == seen.end()) {
877:                   const arrow_type newArrow(*c_iter, arrow.source);
878:                   int              pointO = orientation->restrictPoint(newArrow)[0];

880:                   seen.insert(*c_iter);
881:                   cone->push_back(oriented_arrow_type(newArrow, o));
882:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
883:                 }
884:               }
885:             } else {
886:               const int numSkip = size + o;
887:               int       count   = 0;

889:               for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
890:                 if (count < numSkip) continue;
891:                 if (seen.find(*c_iter) == seen.end()) {
892:                   const arrow_type newArrow(*c_iter, arrow.source);
893:                   int              pointO = orientation->restrictPoint(newArrow)[0];

895:                   seen.insert(*c_iter);
896:                   cone->push_back(oriented_arrow_type(newArrow, o));
897:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
898:                 }
899:               }
900:               count = 0;
901:               for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
902:                 if (count >= numSkip) break;
903:                 if (seen.find(*c_iter) == seen.end()) {
904:                   const arrow_type newArrow(*c_iter, arrow.source);
905:                   int              pointO = orientation->restrictPoint(newArrow)[0];

907:                   seen.insert(*c_iter);
908:                   cone->push_back(oriented_arrow_type(newArrow, o));
909:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
910:                 }
911:               }
912:             }
913:           } else {
914:             if (o == 1) {
915:               for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
916:                 if (seen.find(*c_iter) == seen.end()) {
917:                   const arrow_type newArrow(*c_iter, arrow.source);

919:                   seen.insert(*c_iter);
920:                   cone->push_back(oriented_arrow_type(newArrow, o));
921:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
922:                 }
923:               }
924:             } else {
925:               const int numSkip = o-1;
926:               int       count   = 0;

928:               for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
929:                 if (count < numSkip) continue;
930:                 if (seen.find(*c_iter) == seen.end()) {
931:                   const arrow_type newArrow(*c_iter, arrow.source);

933:                   seen.insert(*c_iter);
934:                   cone->push_back(oriented_arrow_type(newArrow, o));
935:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
936:                 }
937:               }
938:               count = 0;
939:               for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
940:                 if (count >= numSkip) break;
941:                 if (seen.find(*c_iter) == seen.end()) {
942:                   const arrow_type newArrow(*c_iter, arrow.source);

944:                   seen.insert(*c_iter);
945:                   cone->push_back(oriented_arrow_type(newArrow, o));
946:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
947:                 }
948:               }
949:             }
950:           }
951:         }
952:       }
953: #endif
954:     }
955:     template<typename Visitor>
956:     static void orientedStar(const Sieve& sieve, const point_type& p, Visitor& v) {
957:       typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
958:       Retriever sV[2] = {Retriever(200,v), Retriever(200,v)};
959:       int       s     = 0;

961:       v.visitPoint(p, 0);
962:       // Support is guarateed to be ordered correctly
963:       sieve.orientedSupport(p, sV[s]);

965:       while(sV[s].getOrientedSize()) {
966:         const typename Retriever::oriented_point_type *support     = sV[s].getOrientedPoints();
967:         const int                                      supportSize = sV[s].getOrientedSize();
968:         s = 1 - s;

970:         for(int p = 0; p < supportSize; ++p) {
971:           const typename Retriever::point_type& point        = support[p].first;
972:           int                                   pO           = support[p].second == 0 ? 1 : support[p].second;
973:           const int                             pSupportSize = sieve.getSupportSize(point);

975:           if (pO < 0) {
976:             if (pO == -pSupportSize) {
977:               sieve.orientedReverseSupport(point, sV[s]);
978:             } else {
979:               const int numSkip = sieve.getSupportSize(point) + pO;

981:               sV[s].setSkip(sV[s].getSize()+numSkip);
982:               sV[s].setLimit(sV[s].getSize()+pSupportSize);
983:               sieve.orientedReverseSupport(point, sV[s]);
984:               sieve.orientedReverseSupport(point, sV[s]);
985:               sV[s].setSkip(0);
986:               sV[s].setLimit(0);
987:             }
988:           } else {
989:             if (pO == 1) {
990:               sieve.orientedSupport(point, sV[s]);
991:             } else {
992:               const int numSkip = pO-1;

994:               sV[s].setSkip(sV[s].getSize()+numSkip);
995:               sV[s].setLimit(sV[s].getSize()+pSupportSize);
996:               sieve.orientedSupport(point, sV[s]);
997:               sieve.orientedSupport(point, sV[s]);
998:               sV[s].setSkip(0);
999:               sV[s].setLimit(0);
1000:             }
1001:           }
1002:         }
1003:         sV[1-s].clear();
1004:       }
1005:     }
1006:   };

1008:   namespace IFSieveDef {
1009:     template<typename PointType_>
1010:     class Sequence {
1011:     public:
1012:       typedef PointType_ point_type;
1013:       class const_iterator {
1014:       public:
1015:         // Standard iterator typedefs
1016:         typedef std::input_iterator_tag iterator_category;
1017:         typedef PointType_              value_type;
1018:         typedef int                     difference_type;
1019:         typedef value_type*             pointer;
1020:         typedef value_type&             reference;
1021:       protected:
1022:         const point_type *_data;
1023:         int               _pos;
1024:       public:
1025:         const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
1026:         virtual ~const_iterator() {};
1027:       public:
1028:         virtual bool              operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
1029:         virtual bool              operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
1030:         virtual const value_type  operator*() const {return this->_data[this->_pos];};
1031:         virtual const_iterator&   operator++() {++this->_pos; return *this;};
1032:         virtual const_iterator    operator++(int n) {
1033:           const_iterator tmp(*this);
1034:           ++this->_pos;
1035:           return tmp;
1036:         };
1037:       };
1038:       typedef const_iterator iterator;
1039:     protected:
1040:       const point_type *_data;
1041:       int               _begin;
1042:       int               _end;
1043:     public:
1044:       Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
1045:       virtual ~Sequence() {};
1046:     public:
1047:       virtual iterator begin() const {return iterator(_data, _begin);};
1048:       virtual iterator end()   const {return iterator(_data, _end);};
1049:       size_t size() const {return(_end - _begin);}
1050:       void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
1051:     };
1052:   }

1054:   // Interval Final Sieve
1055:   // This is just two CSR matrices that give cones and supports
1056:   //   It is completely static and cannot be resized
1057:   //   It will operator on visitors, rather than sequences (which are messy)
1058:   template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1059:   class IFSieve : public ParallelObject {
1060:   public:
1061:     // Types
1062:     typedef IFSieve<Point_,Allocator_>         this_type;
1063:     typedef Point_                             point_type;
1064:     typedef SimpleArrow<point_type,point_type> arrow_type;
1065:     typedef typename arrow_type::source_type   source_type;
1066:     typedef typename arrow_type::target_type   target_type;
1067:     typedef int                                index_type;
1068:     // Allocators
1069:     typedef Allocator_                                                        point_allocator_type;
1070:     typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1071:     typedef typename point_allocator_type::template rebind<int>::other        int_allocator_type;
1072:     // Interval
1073:     typedef Interval<point_type, point_allocator_type> chart_type;
1074:     // Dynamic structure
1075:     typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1076:     // Iterator interface
1077:     typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1078:     typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1079:     // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1080:     typedef std::set<point_type>   pointSet;
1081:     typedef ALE::array<point_type> pointArray;
1082:     typedef pointSet               coneSet;
1083:     typedef pointSet               supportSet;
1084:     typedef pointArray             coneArray;
1085:     typedef pointArray             supportArray;
1086:   protected:
1087:     // Arrow Containers
1088:     typedef index_type *offsets_type;
1089:     typedef point_type *cones_type;
1090:     typedef point_type *supports_type;
1091:     // Decorators
1092:     typedef int        *orientations_type;
1093:   protected:
1094:     // Data
1095:     bool                 indexAllocated;
1096:     offsets_type         coneOffsets;
1097:     offsets_type         supportOffsets;
1098:     bool                 pointAllocated;
1099:     index_type           maxConeSize;
1100:     index_type           maxSupportSize;
1101:     index_type           baseSize;
1102:     index_type           capSize;
1103:     cones_type           cones;
1104:     supports_type        supports;
1105:     bool                 orientCones;
1106:     orientations_type    coneOrientations;
1107:     chart_type           chart;
1108:     int_allocator_type   intAlloc;
1109:     index_allocator_type indexAlloc;
1110:     point_allocator_type pointAlloc;
1111:     newpoints_type       newCones;
1112:     newpoints_type       newSupports;
1113:     // Sequences
1114:     Obj<coneSequence>    coneSeq;
1115:     Obj<supportSequence> supportSeq;
1116:   protected: // Memory Management
1117:     void createIndices() {
1118:       this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1119:       this->coneOffsets -= this->chart.min();
1120:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1121:       this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1122:       this->supportOffsets -= this->chart.min();
1123:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1124:       this->indexAllocated = true;
1125:     };
1126:     void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1127:       if (*coneOffsets) {
1128:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1129:         *coneOffsets += chart.min();
1130:         indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1131:         *coneOffsets = NULL;
1132:       }
1133:       if (*supportOffsets) {
1134:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1135:         *supportOffsets += chart.min();
1136:         indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1137:         *supportOffsets = NULL;
1138:       }
1139:     };
1140:     void destroyIndices() {
1141:       this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1142:       this->indexAllocated = false;
1143:       this->maxConeSize    = -1;
1144:       this->maxSupportSize = -1;
1145:       this->baseSize       = -1;
1146:       this->capSize        = -1;
1147:     };
1148:     void createPoints() {
1149:       this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1150:       for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1151:       this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1152:       for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1153:       if (orientCones) {
1154:         this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1155:         for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1156:       }
1157:       this->pointAllocated = true;
1158:     };
1159:     void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1160:       if (*cones) {
1161:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1162:         pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1163:         *cones = NULL;
1164:       }
1165:       if (*supports) {
1166:         for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1167:         pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1168:         *supports = NULL;
1169:       }
1170:       if (*coneOrientations) {
1171:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1172:         intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1173:         *coneOrientations = NULL;
1174:       }
1175:     };
1176:     void destroyPoints() {
1177:       this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1178:       this->pointAllocated = false;
1179:     };
1180:     void prefixSum(const offsets_type array) {
1181:       for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1182:         array[p] = array[p] + array[p-1];
1183:       }
1184:     };
1185:     void calculateBaseAndCapSize() {
1186:       this->baseSize = 0;
1187:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1188:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1189:           ++this->baseSize;
1190:         }
1191:       }
1192:       this->capSize = 0;
1193:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1194:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1195:           ++this->capSize;
1196:         }
1197:       }
1198:     };
1199:   public:
1200:     IFSieve(const MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), indexAllocated(false), coneOffsets(NULL), supportOffsets(NULL), pointAllocated(false), maxConeSize(-1), maxSupportSize(-1), baseSize(-1), capSize(-1), cones(NULL), supports(NULL), orientCones(true), coneOrientations(NULL) {
1201:       this->coneSeq    = new coneSequence(NULL, 0, 0);
1202:       this->supportSeq = new supportSequence(NULL, 0, 0);
1203:     };
1204:     IFSieve(const MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : ParallelObject(comm, debug), indexAllocated(false), coneOffsets(NULL), supportOffsets(NULL), pointAllocated(false), maxConeSize(-1), maxSupportSize(-1), baseSize(-1), capSize(-1), cones(NULL), supports(NULL), orientCones(true), coneOrientations(NULL) {
1205:       this->setChart(chart_type(min, max));
1206:       this->coneSeq    = new coneSequence(NULL, 0, 0);
1207:       this->supportSeq = new supportSequence(NULL, 0, 0);
1208:     };
1209:     ~IFSieve() {
1210:       this->destroyPoints();
1211:       this->destroyIndices();
1212:     };
1213:   public: // Accessors
1214:     const chart_type& getChart() const {return this->chart;};
1215:     void setChart(const chart_type& chart) {
1216:       this->destroyPoints();
1217:       this->destroyIndices();
1218:       this->chart = chart;
1219:       this->createIndices();
1220:     };
1221:     index_type getMaxConeSize() const {
1222:       return this->maxConeSize;
1223:     };
1224:     index_type getMaxSupportSize() const {
1225:       return this->maxSupportSize;
1226:     };
1227:     bool baseContains(const point_type& p) const {
1228:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1229:       this->chart.checkPoint(p);

1231:       if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1232:         return true;
1233:       }
1234:       return false;
1235:     };
1236:     bool orientedCones() const {return this->orientCones;};
1237:   public: // Construction
1238:     index_type getConeSize(const point_type& p) const {
1239:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1240:       this->chart.checkPoint(p);
1241:       return this->coneOffsets[p+1]-this->coneOffsets[p];
1242:     };
1243:     void setConeSize(const point_type& p, const index_type c) {
1244:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1245:       this->chart.checkPoint(p);
1246:       this->coneOffsets[p+1] = c;
1247:       this->maxConeSize = std::max(this->maxConeSize, c);
1248:     };
1249:     void addConeSize(const point_type& p, const index_type c) {
1250:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1251:       this->chart.checkPoint(p);
1252:       this->coneOffsets[p+1] += c;
1253:       this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1254:     };
1255:     index_type getSupportSize(const point_type& p) const {
1256:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1257:       this->chart.checkPoint(p);
1258:       return this->supportOffsets[p+1]-this->supportOffsets[p];
1259:     };
1260:     void setSupportSize(const point_type& p, const index_type s) {
1261:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1262:       this->chart.checkPoint(p);
1263:       this->supportOffsets[p+1] = s;
1264:       this->maxSupportSize = std::max(this->maxSupportSize, s);
1265:     };
1266:     void addSupportSize(const point_type& p, const index_type s) {
1267:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1268:       this->chart.checkPoint(p);
1269:       this->supportOffsets[p+1] += s;
1270:       this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1271:     };
1272:     void allocate() {
1273:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1274:       this->prefixSum(this->coneOffsets);
1275:       this->prefixSum(this->supportOffsets);
1276:       this->createPoints();
1277:       this->calculateBaseAndCapSize();
1278:     }
1279:     // Purely for backwards compatibility
1280:     template<typename Color>
1281:     void addArrow(const point_type& p, const point_type& q, const Color c, const bool forceSupport = false) {
1282:       this->addArrow(p, q, forceSupport);
1283:     }
1284:     void addArrow(const point_type& p, const point_type& q, const bool forceSupport = false) {
1285:       if (!this->chart.hasPoint(q)) {
1286:         if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1287:           const index_type start = this->coneOffsets[q];
1288:           const index_type end   = this->coneOffsets[q+1];

1290:           for(int c = start; c < end; ++c) {
1291:             this->newCones[q].push_back(this->cones[c]);
1292:           }
1293:         }
1294:         this->newCones[q].push_back(p);
1295:       }
1296:       if (!this->chart.hasPoint(p) || forceSupport) {
1297:         if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1298:           const index_type start = this->supportOffsets[p];
1299:           const index_type end   = this->supportOffsets[p+1];

1301:           for(int s = start; s < end; ++s) {
1302:             this->newSupports[p].push_back(this->supports[s]);
1303:           }
1304:         }
1305:         this->newSupports[p].push_back(q);
1306:       }
1307:     };
1308:     void reallocate() {
1309:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1310:       if (!this->newCones.size() && !this->newSupports.size()) return;
1311:       const chart_type     oldChart            = this->chart;
1312:       offsets_type         oldConeOffsets      = this->coneOffsets;
1313:       offsets_type         oldSupportOffsets   = this->supportOffsets;
1314:       cones_type           oldCones            = this->cones;
1315:       supports_type        oldSupports         = this->supports;
1316:       orientations_type    oldConeOrientations = this->coneOrientations;
1317:       point_type           min                 = this->chart.min();
1318:       point_type           max                 = this->chart.max()-1;

1320:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1321:         min = std::min(min, c_iter->first);
1322:         max = std::max(max, c_iter->first);
1323:       }
1324:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1325:         min = std::min(min, s_iter->first);
1326:         max = std::max(max, s_iter->first);
1327:       }
1328:       this->chart = chart_type(min, max+1);
1329:       this->createIndices();
1330:       // Copy sizes (converted from offsets)
1331:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1332:         this->coneOffsets[p+1]    = oldConeOffsets[p+1]-oldConeOffsets[p];
1333:         this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1334:       }
1335:       // Inject new sizes
1336:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1337:         this->coneOffsets[c_iter->first+1]    = c_iter->second.size();
1338:         this->maxConeSize                     = std::max(this->maxConeSize,    (int) c_iter->second.size());
1339:       }
1340:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1341:         this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1342:         this->maxSupportSize                  = std::max(this->maxSupportSize, (int) s_iter->second.size());
1343:       }
1344:       this->prefixSum(this->coneOffsets);
1345:       this->prefixSum(this->supportOffsets);
1346:       this->createPoints();
1347:       this->calculateBaseAndCapSize();
1348:       // Copy cones and supports
1349:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1350:         const index_type cStart  = this->coneOffsets[p];
1351:         const index_type cOStart = oldConeOffsets[p];
1352:         const index_type cOEnd   = oldConeOffsets[p+1];
1353:         const index_type sStart  = this->supportOffsets[p];
1354:         const index_type sOStart = oldSupportOffsets[p];
1355:         const index_type sOEnd   = oldSupportOffsets[p+1];

1357:         for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1358:           this->cones[c] = oldCones[cO];
1359:         }
1360:         for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1361:           this->supports[s] = oldSupports[sO];
1362:         }
1363:         if (this->orientCones) {
1364:           for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1365:             this->coneOrientations[c] = oldConeOrientations[cO];
1366:           }
1367:         }
1368:       }
1369:       // Inject new cones and supports
1370:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1371:         index_type start = this->coneOffsets[c_iter->first];

1373:         for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1374:           this->cones[start++] = *p_iter;
1375:         }
1376:         if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1377:       }
1378:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1379:         index_type start = this->supportOffsets[s_iter->first];

1381:         for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1382:           this->supports[start++] = *p_iter;
1383:         }
1384:         if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1385:       }
1386:       this->newCones.clear();
1387:       this->newSupports.clear();
1388:       this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1389:       this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1390:     };
1391:     // Recalculate the support offsets and fill the supports
1392:     //   This is used when an IFSieve is being used as a label
1393:     void recalculateLabel() {
1394:       ISieveVisitor::PointRetriever<IFSieve> v(1);

1396:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1397:         this->supportOffsets[p+1] = 0;
1398:       }
1399:       this->maxSupportSize = 0;
1400:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1401:         this->cone(p, v);
1402:         if (v.getSize()) {
1403:           this->supportOffsets[v.getPoints()[0]+1]++;
1404:           this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1405:         }
1406:         v.clear();
1407:       }
1408:       this->prefixSum(this->supportOffsets);
1409:       this->calculateBaseAndCapSize();
1410:       this->symmetrize();
1411:     };
1412:     void setCone(const point_type cone[], const point_type& p) {
1413:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1414:       this->chart.checkPoint(p);
1415:       const index_type start = this->coneOffsets[p];
1416:       const index_type end   = this->coneOffsets[p+1];

1418:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1419:         this->cones[c] = cone[i];
1420:       }
1421:     };
1422:     void setCone(const point_type cone, const point_type& p) {
1423:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1424:       this->chart.checkPoint(p);
1425:       const index_type start = this->coneOffsets[p];
1426:       const index_type end   = this->coneOffsets[p+1];

1428:       if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1429:       this->cones[start] = cone;
1430:     };
1431: #if 0
1432:     template<typename PointSequence>
1433:     void setCone(const PointSequence& cone, const point_type& p) {
1434:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1435:       this->chart.checkPoint(p);
1436:       const index_type start = this->coneOffsets[p];
1437:       const index_type end   = this->coneOffsets[p+1];
1438:       if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1439:       typename PointSequence::iterator c_iter = cone.begin();

1441:       for(index_type c = start; c < end; ++c, ++c_iter) {
1442:         this->cones[c] = *c_iter;
1443:       }
1444:     };
1445: #endif
1446:     void setSupport(const point_type& p, const point_type support[]) {
1447:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1448:       this->chart.checkPoint(p);
1449:       const index_type start = this->supportOffsets[p];
1450:       const index_type end   = this->supportOffsets[p+1];

1452:       for(index_type s = start, i = 0; s < end; ++s, ++i) {
1453:         this->supports[s] = support[i];
1454:       }
1455:     };
1456: #if 0
1457:     template<typename PointSequence>
1458:     void setSupport(const point_type& p, const PointSequence& support) {
1459:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1460:       this->chart.checkPoint(p);
1461:       const index_type start = this->supportOffsets[p];
1462:       const index_type end   = this->supportOffsets[p+1];
1463:       if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1464:       typename PointSequence::iterator s_iter = support.begin();

1466:       for(index_type s = start; s < end; ++s, ++s_iter) {
1467:         this->supports[s] = *s_iter;
1468:       }
1469:     };
1470: #endif
1471:     void setConeOrientation(const int coneOrientation[], const point_type& p) {
1472:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1473:       this->chart.checkPoint(p);
1474:       const index_type start = this->coneOffsets[p];
1475:       const index_type end   = this->coneOffsets[p+1];

1477:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1478:         this->coneOrientations[c] = coneOrientation[i];
1479:       }
1480:     };
1481:     void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1482:       for(point_type p = 0; p < numCells; ++p) {
1483:         const index_type start = p*numCorners;
1484:         const index_type end   = (p+1)*numCorners;

1486:         for(index_type c = start; c < end; ++c) {
1487:           const point_type q = cones[c]+offset;

1489:           this->supportOffsets[q+1]++;
1490:         }
1491:       }
1492:       for(point_type p = numCells; p < this->chart.max(); ++p) {
1493:         this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1494:       }
1495:     };
1496:     void symmetrize() {
1497:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1498:       offsets -= this->chart.min();
1499:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1500:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1502:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1503:         const index_type start = this->coneOffsets[p];
1504:         const index_type end   = this->coneOffsets[p+1];

1506:         for(index_type c = start; c < end; ++c) {
1507:           const point_type q = this->cones[c];

1509:           this->chart.checkPoint(q);
1510:           this->supports[this->supportOffsets[q]+offsets[q]] = p;
1511:           ++offsets[q];
1512:         }
1513:       }
1514:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1515:       indexAlloc.deallocate(offsets, this->chart.size()+1);
1516:     }
1517:     index_type getBaseSize() const {
1518:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1519:       return this->baseSize;
1520:     }
1521:     index_type getCapSize() const {
1522:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1523:       return this->capSize;
1524:     }
1525:   public: // Traversals
1526:     template<typename Visitor>
1527:     void roots(const Visitor& v) const {
1528:       this->roots(const_cast<Visitor&>(v));
1529:     }
1530:     template<typename Visitor>
1531:     void roots(Visitor& v) const {
1532:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1534:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1535:         if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1536:           if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1537:             v.visitPoint(p);
1538:           }
1539:         }
1540:       }
1541:     }
1542:     template<typename Visitor>
1543:     void leaves(const Visitor& v) const {
1544:       this->leaves(const_cast<Visitor&>(v));
1545:     }
1546:     template<typename Visitor>
1547:     void leaves(Visitor& v) const {
1548:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1550:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1551:         if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1552:           if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1553:             v.visitPoint(p);
1554:           }
1555:         }
1556:       }
1557:     }
1558:     template<typename Visitor>
1559:     void base(const Visitor& v) const {
1560:       this->base(const_cast<Visitor&>(v));
1561:     }
1562:     template<typename Visitor>
1563:     void base(Visitor& v) const {
1564:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1566:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1567:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1568:           v.visitPoint(p);
1569:         }
1570:       }
1571:     }
1572:     template<typename Visitor>
1573:     void cap(const Visitor& v) const {
1574:       this->cap(const_cast<Visitor&>(v));
1575:     }
1576:     template<typename Visitor>
1577:     void cap(Visitor& v) const {
1578:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1580:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1581:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1582:           v.visitPoint(p);
1583:         }
1584:       }
1585:     }
1586:     template<typename PointSequence, typename Visitor>
1587:     void cone(const PointSequence& points, Visitor& v) const {
1588:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1589:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1590:         const point_type p = *p_iter;
1591:         this->chart.checkPoint(p);
1592:         const index_type start = this->coneOffsets[p];
1593:         const index_type end   = this->coneOffsets[p+1];

1595:         for(index_type c = start; c < end; ++c) {
1596:           v.visitArrow(arrow_type(this->cones[c], p));
1597:           v.visitPoint(this->cones[c]);
1598:         }
1599:       }
1600:     }
1601:     template<typename Visitor>
1602:     void cone(const point_type& p, Visitor& v) const {
1603:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1604:       this->chart.checkPoint(p);
1605:       const index_type start = this->coneOffsets[p];
1606:       const index_type end   = this->coneOffsets[p+1];

1608:       for(index_type c = start; c < end; ++c) {
1609:         v.visitArrow(arrow_type(this->cones[c], p));
1610:         v.visitPoint(this->cones[c]);
1611:       }
1612:     }
1613:     const Obj<coneSequence>& cone(const point_type& p) const {
1614:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1615:       if (!this->chart.hasPoint(p)) {
1616:         this->coneSeq->reset(this->cones, 0, 0);
1617:       } else {
1618:         this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1619:       }
1620:       return this->coneSeq;
1621:     }
1622:     template<typename PointSequence, typename Visitor>
1623:     void support(const PointSequence& points, Visitor& v) const {
1624:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1625:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1626:         const point_type p = *p_iter;
1627:         this->chart.checkPoint(p);
1628:         const index_type start = this->supportOffsets[p];
1629:         const index_type end   = this->supportOffsets[p+1];

1631:         for(index_type s = start; s < end; ++s) {
1632:           v.visitArrow(arrow_type(p, this->supports[s]));
1633:           v.visitPoint(this->supports[s]);
1634:         }
1635:       }
1636:     }
1637:     template<typename Visitor>
1638:     void support(const point_type& p, Visitor& v) const {
1639:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1640:       this->chart.checkPoint(p);
1641:       const index_type start = this->supportOffsets[p];
1642:       const index_type end   = this->supportOffsets[p+1];

1644:       for(index_type s = start; s < end; ++s) {
1645:         v.visitArrow(arrow_type(p, this->supports[s]));
1646:         v.visitPoint(this->supports[s]);
1647:       }
1648:     }
1649:     const Obj<supportSequence>& support(const point_type& p) const {
1650:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1651:       if (!this->chart.hasPoint(p)) {
1652:         this->supportSeq->reset(this->supports, 0, 0);
1653:       } else {
1654:         this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1655:       }
1656:       return this->supportSeq;
1657:     }
1658:     template<typename Visitor>
1659:     void orientedCone(const point_type& p, Visitor& v) const {
1660:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1661:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1662:       this->chart.checkPoint(p);
1663:       const index_type start = this->coneOffsets[p];
1664:       const index_type end   = this->coneOffsets[p+1];

1666:       for(index_type c = start; c < end; ++c) {
1667:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1668:         v.visitPoint(this->cones[c], this->coneOrientations[c]);
1669:       }
1670:     }
1671:     template<typename Visitor>
1672:     void orientedReverseCone(const point_type& p, Visitor& v) const {
1673:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1674:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1675:       this->chart.checkPoint(p);
1676:       const index_type start = this->coneOffsets[p];
1677:       const index_type end   = this->coneOffsets[p+1];

1679:       for(index_type c = end-1; c >= start; --c) {
1680:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1681:         v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1682:       }
1683:     }
1684:     template<typename Visitor>
1685:     void orientedSupport(const point_type& p, Visitor& v) const {
1686:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1687:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1688:       this->chart.checkPoint(p);
1689:       const index_type start = this->supportOffsets[p];
1690:       const index_type end   = this->supportOffsets[p+1];

1692:       for(index_type s = start; s < end; ++s) {
1693:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1694:         //v.visitPoint(this->supports[s], this->supportOrientations[s]);
1695:         v.visitArrow(arrow_type(this->supports[s], p), 0);
1696:         v.visitPoint(this->supports[s], 0);
1697:       }
1698:     }
1699:     template<typename Visitor>
1700:     void orientedReverseSupport(const point_type& p, Visitor& v) const {
1701:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1702:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1703:       this->chart.checkPoint(p);
1704:       const index_type start = this->supportOffsets[p];
1705:       const index_type end   = this->supportOffsets[p+1];

1707:       for(index_type s = end-1; s >= start; --s) {
1708:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1709:         //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
1710:         v.visitArrow(arrow_type(this->supports[s], p), 0);
1711:         v.visitPoint(this->supports[s], 0);
1712:       }
1713:     }
1714:     // Currently does only 1 level
1715:     //   Does not check for uniqueness
1716:     template<typename Visitor>
1717:     void meet(const point_type& p, const point_type& q, Visitor& v) const {
1718:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1719:       this->chart.checkPoint(p);
1720:       this->chart.checkPoint(q);
1721:       const index_type startP = this->coneOffsets[p];
1722:       const index_type endP   = this->coneOffsets[p+1];
1723:       const index_type startQ = this->coneOffsets[q];
1724:       const index_type endQ   = this->coneOffsets[q+1];

1726:       for(index_type cP = startP; cP < endP; ++cP) {
1727:         const point_type& c1 = this->cones[cP];

1729:         for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1730:           if (c1 == this->cones[cQ]) v.visitPoint(c1);
1731:         }
1732:         if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1733:       }
1734:     }
1735:     // Currently does only 1 level
1736:     template<typename Sequence, typename Visitor>
1737:     void join(const Sequence& points, Visitor& v) const {
1738:       typedef std::set<point_type> pointSet;
1739:       pointSet intersect[2] = {pointSet(), pointSet()};
1740:       pointSet tmp;
1741:       int      p = 0;
1742:       int      c = 0;

1744:       for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1745:         this->chart.checkPoint(*p_iter);
1746:         tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1747:         if (p == 0) {
1748:           intersect[1-c].insert(tmp.begin(), tmp.end());
1749:           p++;
1750:         } else {
1751:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1752:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1753:           intersect[c].clear();
1754:         }
1755:         c = 1 - c;
1756:         tmp.clear();
1757:       }
1758:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1759:         v.visitPoint(*p_iter);
1760:       }
1761:     }
1762:   public: // Viewing
1763:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1764:       ostringstream txt;
1765:       int rank;

1767:       if (comm == MPI_COMM_NULL) {
1768:         comm = this->comm();
1769:         rank = this->commRank();
1770:       } else {
1771:         MPI_Comm_rank(comm, &rank);
1772:       }
1773:       if (name == "") {
1774:         if(rank == 0) {
1775:           txt << "viewing an IFSieve" << std::endl;
1776:         }
1777:       } else {
1778:         if(rank == 0) {
1779:           txt << "viewing IFSieve '" << name << "'" << std::endl;
1780:         }
1781:       }
1782:       if(rank == 0) {
1783:         txt << "cap --> base:" << std::endl;
1784:       }
1785:       ISieveVisitor::PrintVisitor pV(txt, rank);
1786:       this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1787:       PetscSynchronizedPrintf(comm, txt.str().c_str());
1788:       PetscSynchronizedFlush(comm);
1789:       ostringstream txt2;

1791:       if(rank == 0) {
1792:         txt2 << "base <-- cap:" << std::endl;
1793:       }
1794:       ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1795:       this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1796:       PetscSynchronizedPrintf(comm, txt2.str().c_str());
1797:       PetscSynchronizedFlush(comm);
1798:       if (orientCones) {
1799:         ostringstream txt3;

1801:         if(rank == 0) {
1802:           txt3 << "Orientation:" << std::endl;
1803:         }
1804:         ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1805:         this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1806:         PetscSynchronizedPrintf(comm, txt3.str().c_str());
1807:         PetscSynchronizedFlush(comm);
1808:       }
1809:     };
1810:   };

1812:   class ISieveConverter {
1813:   public:
1814:     template<typename Sieve, typename ISieve, typename Renumbering>
1815:     static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1816:       // First construct a renumbering of the sieve points
1817:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
1818:       const Obj<typename Sieve::capSequence>&  cap  = sieve.cap();
1819:       typename ISieve::point_type              min  = 0;
1820:       typename ISieve::point_type              max  = 0;

1822:       if (renumber) {
1823:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1824:           renumbering[*b_iter] = max++;
1825:         }
1826:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1827:           if (renumbering.find(*c_iter) == renumbering.end()) {
1828:             renumbering[*c_iter] = max++;
1829:           }
1830:         }
1831:       } else {
1832:         if (base->size()) {
1833:           min = *base->begin();
1834:           max = *base->begin();
1835:         } else if (cap->size()) {
1836:           min = *cap->begin();
1837:           max = *cap->begin();
1838:         }
1839:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1840:           min = std::min(min, *b_iter);
1841:           max = std::max(max, *b_iter);
1842:         }
1843:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1844:           min = std::min(min, *c_iter);
1845:           max = std::max(max, *c_iter);
1846:         }
1847:         if (base->size() || cap->size()) {
1848:           ++max;
1849:         }
1850:         for(typename ISieve::point_type p = min; p < max; ++p) {
1851:           renumbering[p] = p;
1852:         }
1853:       }
1854:       // Create the ISieve
1855:       isieve.setChart(typename ISieve::chart_type(min, max));
1856:       // Set cone and support sizes
1857:       size_t maxSize = 0;

1859:       for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1860:         const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);

1862:         isieve.setConeSize(renumbering[*b_iter], cone->size());
1863:         maxSize = std::max(maxSize, cone->size());
1864:       }
1865:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1866:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);

1868:         isieve.setSupportSize(renumbering[*c_iter], support->size());
1869:         maxSize = std::max(maxSize, support->size());
1870:       }
1871:       isieve.allocate();
1872:       // Fill up cones and supports
1873:       typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];

1875:       for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1876:         const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1877:         int i = 0;

1879:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1880:           points[i] = renumbering[*c_iter];
1881:         }
1882:         isieve.setCone(points, renumbering[*b_iter]);
1883:       }
1884:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1885:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1886:         int i = 0;

1888:         for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
1889:           points[i] = renumbering[*s_iter];
1890:         }
1891:         isieve.setSupport(renumbering[*c_iter], points);
1892:       }
1893:       delete [] points;
1894:     }
1895:     template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
1896:     static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
1897:       if (isieve.getMaxConeSize() < 0) return;
1898:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
1899:       int *orientations = new int[isieve.getMaxConeSize()];

1901:       for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1902:         const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1903:         int i = 0;

1905:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1906:           typename ArrowSection::point_type arrow(*c_iter, *b_iter);

1908:           orientations[i] = orientation->restrictPoint(arrow)[0];
1909:         }
1910:         isieve.setConeOrientation(orientations, renumbering[*b_iter]);
1911:       }
1912:       delete [] orientations;
1913:     }
1914:     template<typename Section, typename ISection, typename Renumbering>
1915:     static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
1916:       const typename Section::chart_type& chart = coordinates.getChart();
1917:       typename ISection::point_type       min   = *chart.begin();
1918:       typename ISection::point_type       max   = *chart.begin();

1920:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1921:         min = std::min(min, *p_iter);
1922:         max = std::max(max, *p_iter);
1923:       }
1924:       icoordinates.setChart(typename ISection::chart_type(min, max+1));
1925:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1926:         icoordinates.setFiberDimension(*p_iter, coordinates.getFiberDimension(*p_iter));
1927:       }
1928:       icoordinates.allocatePoint();
1929:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1930:         icoordinates.updatePoint(*p_iter, coordinates.restrictPoint(*p_iter));
1931:       }
1932:     }
1933:     template<typename IMesh, typename Label>
1934:     static void convertLabel(IMesh& imesh, const std::string& name, const Obj<Label>& oldLabel) {
1935:       const Obj<typename IMesh::label_type>&        label = imesh.createLabel(name);
1936:       const typename IMesh::sieve_type::chart_type& chart = imesh.getSieve()->getChart();
1937:       int                                           size  = 0;

1939:       label->setChart(chart);
1940:       for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1941:         const int coneSize = oldLabel->cone(p)->size();

1943:         label->setConeSize(p, coneSize);
1944:         size += coneSize;
1945:       }
1946:       if (size) {label->setSupportSize(0, size);}
1947:       label->allocate();
1948:       for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1949:         const Obj<typename Label::coneSequence>& cone = oldLabel->cone(p);

1951:         if (cone->size()) {
1952:           label->setCone(*cone->begin(), p);
1953:         }
1954:       }
1955:       label->recalculateLabel();
1956:     }
1957:     template<typename Mesh, typename IMesh, typename Renumbering>
1958:     static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
1959:       convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
1960:       imesh.stratify();
1961:       convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
1962:       convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
1963:       if (mesh.hasRealSection("normals")) {
1964:         convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
1965:       }
1966:       const typename Mesh::labels_type& labels = mesh.getLabels();
1967:       std::string heightName("height");
1968:       std::string depthName("depth");

1970:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
1971: #ifdef IMESH_NEW_LABELS
1972:         if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
1973:           convertLabel(imesh, l_iter->first, l_iter->second);
1974:         }
1975: #else
1976:         imesh.setLabel(l_iter->first, l_iter->second);
1977: #endif
1978:       }
1979:     }
1980:   };

1982:   class ISieveSerializer {
1983:   public:
1984:     template<typename ISieve>
1985:     static void writeSieve(const std::string& filename, ISieve& sieve) {
1986:       std::ofstream fs;

1988:       if (sieve.commRank() == 0) {
1989:         fs.open(filename.c_str());
1990:       }
1991:       writeSieve(fs, sieve);
1992:       if (sieve.commRank() == 0) {
1993:         fs.close();
1994:       }
1995:     };
1996:     template<typename ISieve>
1997:     static void writeSieve(std::ofstream& fs, ISieve& sieve) {
1998:       typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
1999:       const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
2000:       typename ISieve::point_type             min   = chart->min();
2001:       typename ISieve::point_type             max   = chart->max();
2002:       PetscInt                               *mins  = new PetscInt[sieve.commSize()];
2003:       PetscInt                               *maxs  = new PetscInt[sieve.commSize()];

2005:       // Write sizes
2006:       if (sieve.commRank() == 0) {
2007:         // Write local
2008:         fs << min <<" "<< max << std::endl;
2009:         for(typename ISieve::point_type p = min; p < max; ++p) {
2010:           fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
2011:         }
2012:         // Receive and write remote
2013:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2014:           PetscInt       minp, maxp;
2015:           PetscInt      *sizes;
2016:           MPI_Status     status;

2019:           MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2020:           MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2021:           PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2022:           MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2023:           fs << minp <<" "<< maxp << std::endl;
2024:           for(PetscInt s = 0; s < maxp - minp; ++s) {
2025:             fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
2026:           }
2027:           PetscFree(sizes);CHKERRXX(ierr);
2028:           mins[pr] = minp;
2029:           maxs[pr] = maxp;
2030:         }
2031:       } else {
2032:         // Send remote
2033:         //PetscInt       min = chart->min();
2034:         //PetscInt       max = chart->max();
2035:         PetscInt       s   = 0;
2036:         PetscInt      *sizes;

2039:         MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2040:         MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2041:         PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
2042:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2043:           sizes[s*2+0] = sieve.getConeSize(p);
2044:           sizes[s*2+1] = sieve.getSupportSize(p);
2045:         }
2046:         MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2047:         PetscFree(sizes);CHKERRXX(ierr);
2048:       }
2049:       // Write covering
2050:       if (sieve.commRank() == 0) {
2051:         // Write local
2052:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);

2054:         for(typename ISieve::point_type p = min; p < max; ++p) {
2055:           sieve.cone(p, pV);
2056:           const typename Visitor::point_type *cone  = pV.getPoints();
2057:           const int                           cSize = pV.getSize();

2059:           if (cSize > 0) {
2060:             for(int c = 0; c < cSize; ++c) {
2061:               if (c) {fs << " ";}
2062:               fs << cone[c];
2063:             }
2064:             fs << std::endl;
2065:           }
2066:           pV.clear();

2068:           sieve.orientedCone(p, pV);
2069:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2070:           const int                                    oSize = pV.getOrientedSize();

2072:           if (oSize > 0) {
2073:             for(int c = 0; c < oSize; ++c) {
2074:               if (c) {fs << " ";}
2075:               fs << oCone[c].second;
2076:             }
2077:             fs << std::endl;
2078:           }
2079:           pV.clear();

2081:           sieve.support(p, pV);
2082:           const typename Visitor::point_type *support = pV.getPoints();
2083:           const int                           sSize   = pV.getSize();

2085:           if (sSize > 0) {
2086:             for(int s = 0; s < sSize; ++s) {
2087:               if (s) {fs << " ";}
2088:               fs << support[s];
2089:             }
2090:             fs << std::endl;
2091:           }
2092:           pV.clear();
2093:         }
2094:         // Receive and write remote
2095:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2096:           PetscInt       size;
2097:           PetscInt      *data;
2098:           PetscInt       off = 0;
2099:           MPI_Status     status;

2102:           MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2103:           PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2104:           MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2105:           for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2106:             PetscInt cSize = data[off++];

2108:             fs << cSize << std::endl;
2109:             if (cSize > 0) {
2110:               for(int c = 0; c < cSize; ++c) {
2111:                 if (c) {fs << " ";}
2112:                 fs << data[off++];
2113:               }
2114:               fs << std::endl;
2115:             }
2116:             PetscInt oSize = data[off++];

2118:             if (oSize > 0) {
2119:               for(int c = 0; c < oSize; ++c) {
2120:                 if (c) {fs << " ";}
2121:                 fs << data[off++];
2122:               }
2123:               fs << std::endl;
2124:             }
2125:             PetscInt sSize = data[off++];

2127:             fs << sSize << std::endl;
2128:             if (sSize > 0) {
2129:               for(int s = 0; s < sSize; ++s) {
2130:                 if (s) {fs << " ";}
2131:                 fs << data[off++];
2132:               }
2133:               fs << std::endl;
2134:             }
2135:           }
2136:           assert(off == size);
2137:           PetscFree(data);CHKERRXX(ierr);
2138:         }
2139:       } else {
2140:         // Send remote
2141:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2142:         PetscInt totalConeSize    = 0;
2143:         PetscInt totalSupportSize = 0;

2145:         for(typename ISieve::point_type p = min; p < max; ++p) {
2146:           totalConeSize    += sieve.getConeSize(p);
2147:           totalSupportSize += sieve.getSupportSize(p);
2148:         }
2149:         PetscInt       size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
2150:         PetscInt       off  = 0;
2151:         PetscInt      *data;

2154:         MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2155:         // There is no nice way to make a generic MPI type here. Really sucky
2156:         PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2157:         for(typename ISieve::point_type p = min; p < max; ++p) {
2158:           sieve.cone(p, pV);
2159:           const typename Visitor::point_type *cone  = pV.getPoints();
2160:           const int                           cSize = pV.getSize();

2162:           data[off++] = cSize;
2163:           for(int c = 0; c < cSize; ++c) {
2164:             data[off++] = cone[c];
2165:           }
2166:           pV.clear();

2168:           sieve.orientedCone(p, pV);
2169:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2170:           const int                                    oSize = pV.getOrientedSize();

2172:           data[off++] = oSize;
2173:           for(int c = 0; c < oSize; ++c) {
2174:             data[off++] = oCone[c].second;
2175:           }
2176:           pV.clear();

2178:           sieve.support(p, pV);
2179:           const typename Visitor::point_type *support = pV.getPoints();
2180:           const int                           sSize   = pV.getSize();

2182:           data[off++] = sSize;
2183:           for(int s = 0; s < sSize; ++s) {
2184:             data[off++] = support[s];
2185:           }
2186:           pV.clear();
2187:         }
2188:         MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2189:         PetscFree(data);CHKERRXX(ierr);
2190:       }
2191:       delete [] mins;
2192:       delete [] maxs;
2193:       // Output renumbering
2194:     };
2195:     template<typename ISieve>
2196:     static void loadSieve(const std::string& filename, ISieve& sieve) {
2197:       std::ifstream fs;

2199:       if (sieve.commRank() == 0) {
2200:         fs.open(filename.c_str());
2201:       }
2202:       loadSieve(fs, sieve);
2203:       if (sieve.commRank() == 0) {
2204:         fs.close();
2205:       }
2206:     };
2207:     template<typename ISieve>
2208:     static void loadSieve(std::ifstream& fs, ISieve& sieve) {
2209:       typename ISieve::point_type min, max;
2210:       PetscInt                   *mins = new PetscInt[sieve.commSize()];
2211:       PetscInt                   *maxs = new PetscInt[sieve.commSize()];
2212:       PetscInt                   *totalConeSizes    = new PetscInt[sieve.commSize()];
2213:       PetscInt                   *totalSupportSizes = new PetscInt[sieve.commSize()];

2215:       // Load sizes
2216:       if (sieve.commRank() == 0) {
2217:         // Load local
2218:         fs >> min;
2219:         fs >> max;
2220:         sieve.setChart(typename ISieve::chart_type(min, max));
2221:         for(typename ISieve::point_type p = min; p < max; ++p) {
2222:           typename ISieve::index_type coneSize, supportSize;

2224:           fs >> coneSize;
2225:           fs >> supportSize;
2226:           sieve.setConeSize(p, coneSize);
2227:           sieve.setSupportSize(p, supportSize);
2228:         }
2229:         // Load and send remote
2230:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2231:           PetscInt       minp, maxp;
2232:           PetscInt      *sizes;

2235:           fs >> minp;
2236:           fs >> maxp;
2237:           MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2238:           MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2239:           PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2240:           totalConeSizes[pr]    = 0;
2241:           totalSupportSizes[pr] = 0;
2242:           for(PetscInt s = 0; s < maxp - minp; ++s) {
2243:             fs >> sizes[s*2+0];
2244:             fs >> sizes[s*2+1];
2245:             totalConeSizes[pr]    += sizes[s*2+0];
2246:             totalSupportSizes[pr] += sizes[s*2+1];
2247:           }
2248:           MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2249:           PetscFree(sizes);CHKERRXX(ierr);
2250:           mins[pr] = minp;
2251:           maxs[pr] = maxp;
2252:         }
2253:       } else {
2254:         // Load remote
2255:         PetscInt       s   = 0;
2256:         PetscInt      *sizes;
2257:         MPI_Status     status;

2260:         MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2261:         MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2262:         sieve.setChart(typename ISieve::chart_type(min, max));
2263:         PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2264:         MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2265:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2266:           sieve.setConeSize(p, sizes[s*2+0]);
2267:           sieve.setSupportSize(p, sizes[s*2+1]);
2268:         }
2269:         PetscFree(sizes);CHKERRXX(ierr);
2270:       }
2271:       sieve.allocate();
2272:       // Load covering
2273:       if (sieve.commRank() == 0) {
2274:         // Load local
2275:         typename ISieve::index_type  maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
2276:         typename ISieve::point_type *points  = new typename ISieve::point_type[maxSize];

2278:         for(typename ISieve::point_type p = min; p < max; ++p) {
2279:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
2280:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);

2282:           if (coneSize > 0) {
2283:             for(int c = 0; c < coneSize; ++c) {
2284:               fs >> points[c];
2285:             }
2286:             sieve.setCone(points, p);
2287:             if (sieve.orientedCones()) {
2288:               for(int c = 0; c < coneSize; ++c) {
2289:                 fs >> points[c];
2290:               }
2291:               sieve.setConeOrientation(points, p);
2292:             }
2293:           }
2294:           if (supportSize > 0) {
2295:             for(int s = 0; s < supportSize; ++s) {
2296:               fs >> points[s];
2297:             }
2298:             sieve.setSupport(p, points);
2299:           }
2300:         }
2301:         delete [] points;
2302:         // Load and send remote
2303:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2304:           PetscInt       size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
2305:           PetscInt       off  = 0;
2306:           PetscInt      *data;

2309:           MPI_Send(&size, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2310:           // There is no nice way to make a generic MPI type here. Really sucky
2311:           PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2312:           for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2313:             PetscInt coneSize, supportSize;

2315:             fs >> coneSize;
2316:             data[off++] = coneSize;
2317:             if (coneSize > 0) {
2318:               for(int c = 0; c < coneSize; ++c) {
2319:                 fs >> data[off++];
2320:               }
2321:               for(int c = 0; c < coneSize; ++c) {
2322:                 fs >> data[off++];
2323:               }
2324:             }
2325:             fs >> supportSize;
2326:             data[off++] = supportSize;
2327:             if (supportSize > 0) {
2328:               for(int s = 0; s < supportSize; ++s) {
2329:                 fs >> data[off++];
2330:               }
2331:             }
2332:           }
2333:           assert(off == size);
2334:           MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2335:           PetscFree(data);CHKERRXX(ierr);
2336:         }
2337:         delete [] mins;
2338:         delete [] maxs;
2339:       } else {
2340:         // Load remote
2341:         PetscInt       size;
2342:         PetscInt      *data;
2343:         PetscInt       off = 0;
2344:         MPI_Status     status;

2347:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2348:         PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2349:         MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2350:         for(typename ISieve::point_type p = min; p < max; ++p) {
2351:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
2352:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2353:           PetscInt cs = data[off++];

2355:           assert(cs == coneSize);
2356:           if (coneSize > 0) {
2357:             sieve.setCone(&data[off], p);
2358:             off += coneSize;
2359:             if (sieve.orientedCones()) {
2360:               sieve.setConeOrientation(&data[off], p);
2361:               off += coneSize;
2362:             }
2363:           }
2364:           PetscInt ss = data[off++];

2366:           assert(ss == supportSize);
2367:           if (supportSize > 0) {
2368:             sieve.setSupport(p, &data[off]);
2369:             off += supportSize;
2370:           }
2371:         }
2372:         assert(off == size);
2373:       }
2374:       // Load renumbering
2375:     };
2376:   };
2377: }

2379: #endif