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: //#define IMESH_NEW_LABELS

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

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

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

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

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

637:           for(int j = start; j < end; ++j, ++i) {
638:             this->values[i] = j;
639:           }
640:         } else if (!section.getNumSpaces()) {
641:           const int start = this->order.getIndex(p);
642:           const int end   = start + section.getFiberDimension(p);

644:           for(int j = end-1; j >= start; --j, ++i) {
645:             this->values[i] = j;
646:           }
647:         } else {
648:           const int numSpaces = section.getNumSpaces();
649:           int       start     = this->order.getIndex(p);

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

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

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

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

685:           for(int space = 0; space < section.getNumSpaces(); ++space) {
686:             const int  dim = this->section.getFiberDimension(p, space);
687:             const int tDim = this->section.getConstrainedFiberDimension(p, space);
688:             int       cInd = (dim - tDim)-1;

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

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

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

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

802:       while(cV[c].getOrientedSize()) {
803:         const typename Retriever::oriented_point_type *cone     = cV[c].getOrientedPoints();
804:         const int                                      coneSize = cV[c].getOrientedSize();
805:         c = 1 - c;

807:         for(int p = 0; p < coneSize; ++p) {
808:           const typename Retriever::point_type& point     = cone[p].first;
809:           int                                   pO        = cone[p].second == 0 ? 1 : cone[p].second;
810:           const int                             pConeSize = sieve.getConeSize(point);

812:           if (pO < 0) {
813:             if (pO == -pConeSize) {
814:               sieve.orientedReverseCone(point, cV[c]);
815:             } else {
816:               const int numSkip = sieve.getConeSize(point) + pO;

818:               cV[c].setSkip(cV[c].getSize()+numSkip);
819:               cV[c].setLimit(cV[c].getSize()+pConeSize);
820:               sieve.orientedReverseCone(point, cV[c]);
821:               sieve.orientedReverseCone(point, cV[c]);
822:               cV[c].setSkip(0);
823:               cV[c].setLimit(0);
824:             }
825:           } else {
826:             if (pO == 1) {
827:               sieve.orientedCone(point, cV[c]);
828:             } else {
829:               const int numSkip = pO-1;

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

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

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

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

863:           if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
864:             o = -(o+1);
865:           }
866:           if (o < 0) {
867:             const int size = pCone->size();

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

875:                   seen.insert(*c_iter);
876:                   cone->push_back(oriented_arrow_type(newArrow, o));
877:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
878:                 }
879:               }
880:             } else {
881:               const int numSkip = size + o;
882:               int       count   = 0;

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

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

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

914:                   seen.insert(*c_iter);
915:                   cone->push_back(oriented_arrow_type(newArrow, o));
916:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
917:                 }
918:               }
919:             } else {
920:               const int numSkip = o-1;
921:               int       count   = 0;

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

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

939:                   seen.insert(*c_iter);
940:                   cone->push_back(oriented_arrow_type(newArrow, o));
941:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
942:                 }
943:               }
944:             }
945:           }
946:         }
947:       }
948: #endif
949:     };
950:   };

952:   namespace IFSieveDef {
953:     template<typename PointType_>
954:     class Sequence {
955:     public:
956:       typedef PointType_ point_type;
957:       class const_iterator {
958:       public:
959:         // Standard iterator typedefs
960:         typedef std::input_iterator_tag iterator_category;
961:         typedef PointType_              value_type;
962:         typedef int                     difference_type;
963:         typedef value_type*             pointer;
964:         typedef value_type&             reference;
965:       protected:
966:         const point_type *_data;
967:         int               _pos;
968:       public:
969:         const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
970:         ~const_iterator() {};
971:       public:
972:         virtual bool              operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
973:         virtual bool              operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
974:         virtual const value_type  operator*() const {return this->_data[this->_pos];};
975:         virtual const_iterator&   operator++() {++this->_pos; return *this;};
976:         virtual const_iterator    operator++(int n) {
977:           const_iterator tmp(*this);
978:           ++this->_pos;
979:           return tmp;
980:         };
981:       };
982:       typedef const_iterator iterator;
983:     protected:
984:       const point_type *_data;
985:       int               _begin;
986:       int               _end;
987:     public:
988:       Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
989:       ~Sequence() {};
990:     public:
991:       virtual iterator begin() const {return iterator(_data, _begin);};
992:       virtual iterator end()   const {return iterator(_data, _end);};
993:       size_t size() const {return(_end - _begin);}
994:       void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
995:     };
996:   }

998:   // Interval Final Sieve
999:   // This is just two CSR matrices that give cones and supports
1000:   //   It is completely static and cannot be resized
1001:   //   It will operator on visitors, rather than sequences (which are messy)
1002:   template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1003:   class IFSieve : public ParallelObject {
1004:   public:
1005:     // Types
1006:     typedef IFSieve<Point_,Allocator_>         this_type;
1007:     typedef Point_                             point_type;
1008:     typedef SimpleArrow<point_type,point_type> arrow_type;
1009:     typedef typename arrow_type::source_type   source_type;
1010:     typedef typename arrow_type::target_type   target_type;
1011:     typedef int                                index_type;
1012:     // Allocators
1013:     typedef Allocator_                                                        point_allocator_type;
1014:     typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1015:     typedef typename point_allocator_type::template rebind<int>::other        int_allocator_type;
1016:     // Interval
1017:     typedef Interval<point_type, point_allocator_type> chart_type;
1018:     // Dynamic structure
1019:     typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1020:     // Iterator interface
1021:     typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1022:     typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1023:     // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1024:     typedef std::set<point_type>   pointSet;
1025:     typedef ALE::array<point_type> pointArray;
1026:     typedef pointSet               coneSet;
1027:     typedef pointSet               supportSet;
1028:     typedef pointArray             coneArray;
1029:     typedef pointArray             supportArray;
1030:   protected:
1031:     // Arrow Containers
1032:     typedef index_type *offsets_type;
1033:     typedef point_type *cones_type;
1034:     typedef point_type *supports_type;
1035:     // Decorators
1036:     typedef int        *orientations_type;
1037:   protected:
1038:     // Data
1039:     bool                 indexAllocated;
1040:     offsets_type         coneOffsets;
1041:     offsets_type         supportOffsets;
1042:     bool                 pointAllocated;
1043:     index_type           maxConeSize;
1044:     index_type           maxSupportSize;
1045:     index_type           baseSize;
1046:     index_type           capSize;
1047:     cones_type           cones;
1048:     supports_type        supports;
1049:     bool                 orientCones;
1050:     orientations_type    coneOrientations;
1051:     chart_type           chart;
1052:     int_allocator_type   intAlloc;
1053:     index_allocator_type indexAlloc;
1054:     point_allocator_type pointAlloc;
1055:     newpoints_type       newCones;
1056:     newpoints_type       newSupports;
1057:     // Sequences
1058:     Obj<coneSequence>    coneSeq;
1059:     Obj<supportSequence> supportSeq;
1060:   protected: // Memory Management
1061:     void createIndices() {
1062:       this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1063:       this->coneOffsets -= this->chart.min();
1064:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1065:       this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1066:       this->supportOffsets -= this->chart.min();
1067:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1068:       this->indexAllocated = true;
1069:     };
1070:     void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1071:       if (*coneOffsets) {
1072:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1073:         *coneOffsets += chart.min();
1074:         indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1075:         *coneOffsets = NULL;
1076:       }
1077:       if (*supportOffsets) {
1078:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1079:         *supportOffsets += chart.min();
1080:         indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1081:         *supportOffsets = NULL;
1082:       }
1083:     };
1084:     void destroyIndices() {
1085:       this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1086:       this->indexAllocated = false;
1087:       this->maxConeSize    = -1;
1088:       this->maxSupportSize = -1;
1089:       this->baseSize       = -1;
1090:       this->capSize        = -1;
1091:     };
1092:     void createPoints() {
1093:       this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1094:       for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1095:       this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1096:       for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1097:       if (orientCones) {
1098:         this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1099:         for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1100:       }
1101:       this->pointAllocated = true;
1102:     };
1103:     void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1104:       if (*cones) {
1105:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1106:         pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1107:         *cones = NULL;
1108:       }
1109:       if (*supports) {
1110:         for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1111:         pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1112:         *supports = NULL;
1113:       }
1114:       if (*coneOrientations) {
1115:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1116:         intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1117:         *coneOrientations = NULL;
1118:       }
1119:     };
1120:     void destroyPoints() {
1121:       this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1122:       this->pointAllocated = false;
1123:     };
1124:     void prefixSum(const offsets_type array) {
1125:       for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1126:         array[p] = array[p] + array[p-1];
1127:       }
1128:     };
1129:     void calculateBaseAndCapSize() {
1130:       this->baseSize = 0;
1131:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1132:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1133:           ++this->baseSize;
1134:         }
1135:       }
1136:       this->capSize = 0;
1137:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1138:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1139:           ++this->capSize;
1140:         }
1141:       }
1142:     };
1143:   public:
1144:     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) {
1145:       this->coneSeq    = new coneSequence(NULL, 0, 0);
1146:       this->supportSeq = new supportSequence(NULL, 0, 0);
1147:     };
1148:     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) {
1149:       this->setChart(chart_type(min, max));
1150:       this->coneSeq    = new coneSequence(NULL, 0, 0);
1151:       this->supportSeq = new supportSequence(NULL, 0, 0);
1152:     };
1153:     ~IFSieve() {
1154:       this->destroyPoints();
1155:       this->destroyIndices();
1156:     };
1157:   public: // Accessors
1158:     const chart_type& getChart() const {return this->chart;};
1159:     void setChart(const chart_type& chart) {
1160:       this->destroyPoints();
1161:       this->destroyIndices();
1162:       this->chart = chart;
1163:       this->createIndices();
1164:     };
1165:     index_type getMaxConeSize() const {
1166:       return this->maxConeSize;
1167:     };
1168:     index_type getMaxSupportSize() const {
1169:       return this->maxSupportSize;
1170:     };
1171:     bool baseContains(const point_type& p) const {
1172:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1173:       this->chart.checkPoint(p);

1175:       if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1176:         return true;
1177:       }
1178:       return false;
1179:     };
1180:   public: // Construction
1181:     index_type getConeSize(const point_type& p) const {
1182:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1183:       this->chart.checkPoint(p);
1184:       return this->coneOffsets[p+1]-this->coneOffsets[p];
1185:     };
1186:     void setConeSize(const point_type& p, const index_type c) {
1187:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1188:       this->chart.checkPoint(p);
1189:       this->coneOffsets[p+1] = c;
1190:       this->maxConeSize = std::max(this->maxConeSize, c);
1191:     };
1192:     void addConeSize(const point_type& p, const index_type c) {
1193:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1194:       this->chart.checkPoint(p);
1195:       this->coneOffsets[p+1] += c;
1196:       this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1197:     };
1198:     index_type getSupportSize(const point_type& p) const {
1199:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1200:       this->chart.checkPoint(p);
1201:       return this->supportOffsets[p+1]-this->supportOffsets[p];
1202:     };
1203:     void setSupportSize(const point_type& p, const index_type s) {
1204:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1205:       this->chart.checkPoint(p);
1206:       this->supportOffsets[p+1] = s;
1207:       this->maxSupportSize = std::max(this->maxSupportSize, s);
1208:     };
1209:     void addSupportSize(const point_type& p, const index_type s) {
1210:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1211:       this->chart.checkPoint(p);
1212:       this->supportOffsets[p+1] += s;
1213:       this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1214:     };
1215:     void allocate() {
1216:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1217:       this->prefixSum(this->coneOffsets);
1218:       this->prefixSum(this->supportOffsets);
1219:       this->createPoints();
1220:       this->calculateBaseAndCapSize();
1221:     };
1222:     // Purely for backwards compatibility
1223:     template<typename Color>
1224:     void addArrow(const point_type& p, const point_type& q, const Color c) {
1225:       this->addArrow(p, q);
1226:     };
1227:     void addArrow(const point_type& p, const point_type& q) {
1228:       if (!this->chart.hasPoint(q)) {
1229:         if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1230:           const index_type start = this->coneOffsets[q];
1231:           const index_type end   = this->coneOffsets[q+1];

1233:           for(int c = start; c < end; ++c) {
1234:             this->newCones[q].push_back(this->cones[c]);
1235:           }
1236:         }
1237:         this->newCones[q].push_back(p);
1238:       }
1239:       if (!this->chart.hasPoint(p)) {
1240:         if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1241:           const index_type start = this->supportOffsets[p];
1242:           const index_type end   = this->supportOffsets[p+1];

1244:           for(int s = start; s < end; ++s) {
1245:             this->newSupports[p].push_back(this->supports[s]);
1246:           }
1247:         }
1248:         this->newSupports[p].push_back(q);
1249:       }
1250:     };
1251:     void reallocate() {
1252:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1253:       if (!this->newCones.size() && !this->newSupports.size()) return;
1254:       const chart_type     oldChart            = this->chart;
1255:       offsets_type         oldConeOffsets      = this->coneOffsets;
1256:       offsets_type         oldSupportOffsets   = this->supportOffsets;
1257:       cones_type           oldCones            = this->cones;
1258:       supports_type        oldSupports         = this->supports;
1259:       orientations_type    oldConeOrientations = this->coneOrientations;
1260:       point_type           min                 = this->chart.min();
1261:       point_type           max                 = this->chart.max()-1;

1263:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1264:         min = std::min(min, c_iter->first);
1265:         max = std::max(max, c_iter->first);
1266:       }
1267:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1268:         min = std::min(min, s_iter->first);
1269:         max = std::max(max, s_iter->first);
1270:       }
1271:       this->chart = chart_type(min, max+1);
1272:       this->createIndices();
1273:       // Copy sizes (converted from offsets)
1274:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1275:         this->coneOffsets[p+1]    = oldConeOffsets[p+1]-oldConeOffsets[p];
1276:         this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1277:       }
1278:       // Inject new sizes
1279:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1280:         this->coneOffsets[c_iter->first+1]    = c_iter->second.size();
1281:         this->maxConeSize                     = std::max(this->maxConeSize,    (int) c_iter->second.size());
1282:       }
1283:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1284:         this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1285:         this->maxSupportSize                  = std::max(this->maxSupportSize, (int) s_iter->second.size());
1286:       }
1287:       this->prefixSum(this->coneOffsets);
1288:       this->prefixSum(this->supportOffsets);
1289:       this->createPoints();
1290:       this->calculateBaseAndCapSize();
1291:       // Copy cones and supports
1292:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1293:         const index_type cStart  = this->coneOffsets[p];
1294:         const index_type cOStart = oldConeOffsets[p];
1295:         const index_type cOEnd   = oldConeOffsets[p+1];
1296:         const index_type sStart  = this->supportOffsets[p];
1297:         const index_type sOStart = oldSupportOffsets[p];
1298:         const index_type sOEnd   = oldSupportOffsets[p+1];

1300:         for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1301:           this->cones[c] = oldCones[cO];
1302:         }
1303:         for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1304:           this->supports[s] = oldSupports[sO];
1305:         }
1306:         if (this->orientCones) {
1307:           for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1308:             this->coneOrientations[c] = oldConeOrientations[cO];
1309:           }
1310:         }
1311:       }
1312:       // Inject new cones and supports
1313:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1314:         index_type start = this->coneOffsets[c_iter->first];

1316:         for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1317:           this->cones[start++] = *p_iter;
1318:         }
1319:         if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1320:       }
1321:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1322:         index_type start = this->supportOffsets[s_iter->first];

1324:         for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1325:           this->supports[start++] = *p_iter;
1326:         }
1327:         if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1328:       }
1329:       this->newCones.clear();
1330:       this->newSupports.clear();
1331:       this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1332:       this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1333:     };
1334:     // Recalculate the support offsets and fill the supports
1335:     //   This is used when an IFSieve is being used as a label
1336:     void recalculateLabel() {
1337:       ISieveVisitor::PointRetriever<IFSieve> v(1);

1339:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1340:         this->supportOffsets[p+1] = 0;
1341:       }
1342:       this->maxSupportSize = 0;
1343:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1344:         this->cone(p, v);
1345:         if (v.getSize()) {
1346:           this->supportOffsets[v.getPoints()[0]+1]++;
1347:           this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1348:         }
1349:         v.clear();
1350:       }
1351:       this->prefixSum(this->supportOffsets);
1352:       this->calculateBaseAndCapSize();
1353:       this->symmetrize();
1354:     };
1355:     void setCone(const point_type cone[], const point_type& p) {
1356:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1357:       this->chart.checkPoint(p);
1358:       const index_type start = this->coneOffsets[p];
1359:       const index_type end   = this->coneOffsets[p+1];

1361:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1362:         this->cones[c] = cone[i];
1363:       }
1364:     };
1365:     void setCone(const point_type cone, const point_type& p) {
1366:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1367:       this->chart.checkPoint(p);
1368:       const index_type start = this->coneOffsets[p];
1369:       const index_type end   = this->coneOffsets[p+1];

1371:       if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1372:       this->cones[start] = cone;
1373:     };
1374: #if 0
1375:     template<typename PointSequence>
1376:     void setCone(const PointSequence& cone, const point_type& p) {
1377:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1378:       this->chart.checkPoint(p);
1379:       const index_type start = this->coneOffsets[p];
1380:       const index_type end   = this->coneOffsets[p+1];
1381:       if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1382:       typename PointSequence::iterator c_iter = cone.begin();

1384:       for(index_type c = start; c < end; ++c, ++c_iter) {
1385:         this->cones[c] = *c_iter;
1386:       }
1387:     };
1388: #endif
1389:     void setSupport(const point_type& p, const point_type support[]) {
1390:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1391:       this->chart.checkPoint(p);
1392:       const index_type start = this->supportOffsets[p];
1393:       const index_type end   = this->supportOffsets[p+1];

1395:       for(index_type s = start, i = 0; s < end; ++s, ++i) {
1396:         this->supports[s] = support[i];
1397:       }
1398:     };
1399: #if 0
1400:     template<typename PointSequence>
1401:     void setSupport(const point_type& p, const PointSequence& support) {
1402:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1403:       this->chart.checkPoint(p);
1404:       const index_type start = this->supportOffsets[p];
1405:       const index_type end   = this->supportOffsets[p+1];
1406:       if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1407:       typename PointSequence::iterator s_iter = support.begin();

1409:       for(index_type s = start; s < end; ++s, ++s_iter) {
1410:         this->supports[s] = *s_iter;
1411:       }
1412:     };
1413: #endif
1414:     void setConeOrientation(const int coneOrientation[], const point_type& p) {
1415:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1416:       this->chart.checkPoint(p);
1417:       const index_type start = this->coneOffsets[p];
1418:       const index_type end   = this->coneOffsets[p+1];

1420:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1421:         this->coneOrientations[c] = coneOrientation[i];
1422:       }
1423:     };
1424:     void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1425:       for(point_type p = 0; p < numCells; ++p) {
1426:         const index_type start = p*numCorners;
1427:         const index_type end   = (p+1)*numCorners;

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

1432:           this->supportOffsets[q+1]++;
1433:         }
1434:       }
1435:       for(point_type p = numCells; p < this->chart.max(); ++p) {
1436:         this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1437:       }
1438:     };
1439:     void symmetrize() {
1440:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1441:       offsets -= this->chart.min();
1442:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1443:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

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

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

1452:           this->chart.checkPoint(q);
1453:           this->supports[this->supportOffsets[q]+offsets[q]] = p;
1454:           ++offsets[q];
1455:         }
1456:       }
1457:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1458:       indexAlloc.deallocate(offsets, this->chart.size()+1);
1459:     };
1460:     index_type getBaseSize() const {
1461:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1462:       return this->baseSize;
1463:     };
1464:     index_type getCapSize() const {
1465:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1466:       return this->capSize;
1467:     };
1468:   public: // Traversals
1469:     template<typename Visitor>
1470:     void roots(const Visitor& v) const {
1471:       this->roots(const_cast<Visitor&>(v));
1472:     };
1473:     template<typename Visitor>
1474:     void roots(Visitor& v) const {
1475:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1477:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1478:         if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1479:           if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1480:             v.visitPoint(p);
1481:           }
1482:         }
1483:       }
1484:     };
1485:     template<typename Visitor>
1486:     void leaves(const Visitor& v) const {
1487:       this->leaves(const_cast<Visitor&>(v));
1488:     };
1489:     template<typename Visitor>
1490:     void leaves(Visitor& v) const {
1491:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1493:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1494:         if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1495:           if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1496:             v.visitPoint(p);
1497:           }
1498:         }
1499:       }
1500:     };
1501:     template<typename Visitor>
1502:     void base(const Visitor& v) const {
1503:       this->base(const_cast<Visitor&>(v));
1504:     };
1505:     template<typename Visitor>
1506:     void base(Visitor& v) const {
1507:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1509:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1510:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1511:           v.visitPoint(p);
1512:         }
1513:       }
1514:     };
1515:     template<typename Visitor>
1516:     void cap(const Visitor& v) const {
1517:       this->cap(const_cast<Visitor&>(v));
1518:     };
1519:     template<typename Visitor>
1520:     void cap(Visitor& v) const {
1521:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1523:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1524:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1525:           v.visitPoint(p);
1526:         }
1527:       }
1528:     };
1529:     template<typename PointSequence, typename Visitor>
1530:     void cone(const PointSequence& points, Visitor& v) const {
1531:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1532:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1533:         const point_type p = *p_iter;
1534:         this->chart.checkPoint(p);
1535:         const index_type start = this->coneOffsets[p];
1536:         const index_type end   = this->coneOffsets[p+1];

1538:         for(index_type c = start; c < end; ++c) {
1539:           v.visitArrow(arrow_type(this->cones[c], p));
1540:           v.visitPoint(this->cones[c]);
1541:         }
1542:       }
1543:     };
1544:     template<typename Visitor>
1545:     void cone(const point_type& p, Visitor& v) const {
1546:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1547:       this->chart.checkPoint(p);
1548:       const index_type start = this->coneOffsets[p];
1549:       const index_type end   = this->coneOffsets[p+1];

1551:       for(index_type c = start; c < end; ++c) {
1552:         v.visitArrow(arrow_type(this->cones[c], p));
1553:         v.visitPoint(this->cones[c]);
1554:       }
1555:     };
1556:     const Obj<coneSequence>& cone(const point_type& p) const {
1557:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1558:       if (!this->chart.hasPoint(p)) {
1559:         this->coneSeq->reset(this->cones, 0, 0);
1560:       } else {
1561:         this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1562:       }
1563:       return this->coneSeq;
1564:     };
1565:     template<typename PointSequence, typename Visitor>
1566:     void support(const PointSequence& points, Visitor& v) const {
1567:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1568:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1569:         const point_type p = *p_iter;
1570:         this->chart.checkPoint(p);
1571:         const index_type start = this->supportOffsets[p];
1572:         const index_type end   = this->supportOffsets[p+1];

1574:         for(index_type s = start; s < end; ++s) {
1575:           v.visitArrow(arrow_type(p, this->supports[s]));
1576:           v.visitPoint(this->supports[s]);
1577:         }
1578:       }
1579:     };
1580:     template<typename Visitor>
1581:     void support(const point_type& p, Visitor& v) const {
1582:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1583:       this->chart.checkPoint(p);
1584:       const index_type start = this->supportOffsets[p];
1585:       const index_type end   = this->supportOffsets[p+1];

1587:       for(index_type s = start; s < end; ++s) {
1588:         v.visitArrow(arrow_type(p, this->supports[s]));
1589:         v.visitPoint(this->supports[s]);
1590:       }
1591:     };
1592:     const Obj<supportSequence>& support(const point_type& p) const {
1593:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1594:       if (!this->chart.hasPoint(p)) {
1595:         this->supportSeq->reset(this->supports, 0, 0);
1596:       } else {
1597:         this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1598:       }
1599:       return this->supportSeq;
1600:     };
1601:     template<typename Visitor>
1602:     void orientedCone(const point_type& p, Visitor& v) const {
1603:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1604:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1605:       this->chart.checkPoint(p);
1606:       const index_type start = this->coneOffsets[p];
1607:       const index_type end   = this->coneOffsets[p+1];

1609:       for(index_type c = start; c < end; ++c) {
1610:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1611:         v.visitPoint(this->cones[c], this->coneOrientations[c]);
1612:       }
1613:     };
1614:     template<typename Visitor>
1615:     void orientedReverseCone(const point_type& p, Visitor& v) const {
1616:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1617:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1618:       this->chart.checkPoint(p);
1619:       const index_type start = this->coneOffsets[p];
1620:       const index_type end   = this->coneOffsets[p+1];

1622:       for(index_type c = end-1; c >= start; --c) {
1623:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1624:         v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1625:       }
1626:     };
1627:     // Currently does only 1 level
1628:     //   Does not check for uniqueness
1629:     template<typename Visitor>
1630:     void meet(const point_type& p, const point_type& q, Visitor& v) const {
1631:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1632:       this->chart.checkPoint(p);
1633:       this->chart.checkPoint(q);
1634:       const index_type startP = this->coneOffsets[p];
1635:       const index_type endP   = this->coneOffsets[p+1];
1636:       const index_type startQ = this->coneOffsets[q];
1637:       const index_type endQ   = this->coneOffsets[q+1];

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

1642:         for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1643:           if (c1 == this->cones[cQ]) v.visitPoint(c1);
1644:         }
1645:         if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1646:       }
1647:     };
1648:     // Currently does only 1 level
1649:     template<typename Sequence, typename Visitor>
1650:     void join(const Sequence& points, Visitor& v) const {
1651:       typedef std::set<point_type> pointSet;
1652:       pointSet intersect[2] = {pointSet(), pointSet()};
1653:       pointSet tmp;
1654:       int      p = 0;
1655:       int      c = 0;

1657:       for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1658:         this->chart.checkPoint(*p_iter);
1659:         tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1660:         if (p == 0) {
1661:           intersect[1-c].insert(tmp.begin(), tmp.end());
1662:           p++;
1663:         } else {
1664:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1665:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1666:           intersect[c].clear();
1667:         }
1668:         c = 1 - c;
1669:         tmp.clear();
1670:       }
1671:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1672:         v.visitPoint(*p_iter);
1673:       }
1674:     };
1675:   public: // Viewing
1676:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1677:       ostringstream txt;
1678:       int rank;

1680:       if (comm == MPI_COMM_NULL) {
1681:         comm = this->comm();
1682:         rank = this->commRank();
1683:       } else {
1684:         MPI_Comm_rank(comm, &rank);
1685:       }
1686:       if (name == "") {
1687:         if(rank == 0) {
1688:           txt << "viewing an IFSieve" << std::endl;
1689:         }
1690:       } else {
1691:         if(rank == 0) {
1692:           txt << "viewing IFSieve '" << name << "'" << std::endl;
1693:         }
1694:       }
1695:       if(rank == 0) {
1696:         txt << "cap --> base:" << std::endl;
1697:       }
1698:       ISieveVisitor::PrintVisitor pV(txt, rank);
1699:       this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1700:       PetscSynchronizedPrintf(comm, txt.str().c_str());
1701:       PetscSynchronizedFlush(comm);
1702:       ostringstream txt2;

1704:       if(rank == 0) {
1705:         txt2 << "base <-- cap:" << std::endl;
1706:       }
1707:       ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1708:       this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1709:       PetscSynchronizedPrintf(comm, txt2.str().c_str());
1710:       PetscSynchronizedFlush(comm);
1711:       if (orientCones) {
1712:         ostringstream txt3;

1714:         if(rank == 0) {
1715:           txt3 << "Orientation:" << std::endl;
1716:         }
1717:         ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1718:         this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1719:         PetscSynchronizedPrintf(comm, txt3.str().c_str());
1720:         PetscSynchronizedFlush(comm);
1721:       }
1722:     };
1723:   };

1725:   class ISieveConverter {
1726:   public:
1727:     template<typename Sieve, typename ISieve, typename Renumbering>
1728:     static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1729:       // First construct a renumbering of the sieve points
1730:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
1731:       const Obj<typename Sieve::capSequence>&  cap  = sieve.cap();
1732:       typename ISieve::point_type              min  = 0;
1733:       typename ISieve::point_type              max  = 0;

1735:       if (renumber) {
1736:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1737:           renumbering[*b_iter] = max++;
1738:         }
1739:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1740:           if (renumbering.find(*c_iter) == renumbering.end()) {
1741:             renumbering[*c_iter] = max++;
1742:           }
1743:         }
1744:       } else {
1745:         if (base->size()) {
1746:           min = *base->begin();
1747:           max = *base->begin();
1748:         } else if (cap->size()) {
1749:           min = *cap->begin();
1750:           max = *cap->begin();
1751:         }
1752:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1753:           min = std::min(min, *b_iter);
1754:           max = std::max(max, *b_iter);
1755:         }
1756:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1757:           min = std::min(min, *c_iter);
1758:           max = std::max(max, *c_iter);
1759:         }
1760:         if (base->size() || cap->size()) {
1761:           ++max;
1762:         }
1763:         for(typename ISieve::point_type p = min; p < max; ++p) {
1764:           renumbering[p] = p;
1765:         }
1766:       }
1767:       // Create the ISieve
1768:       isieve.setChart(typename ISieve::chart_type(min, max));
1769:       // Set cone and support sizes
1770:       size_t maxSize = 0;

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

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

1781:         isieve.setSupportSize(renumbering[*c_iter], support->size());
1782:         maxSize = std::max(maxSize, support->size());
1783:       }
1784:       isieve.allocate();
1785:       // Fill up cones and supports
1786:       typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];

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

1792:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1793:           points[i] = renumbering[*c_iter];
1794:         }
1795:         isieve.setCone(points, renumbering[*b_iter]);
1796:       }
1797:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1798:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1799:         int i = 0;

1801:         for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
1802:           points[i] = renumbering[*s_iter];
1803:         }
1804:         isieve.setSupport(renumbering[*c_iter], points);
1805:       }
1806:       delete [] points;
1807:     };
1808:     template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
1809:     static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
1810:       if (isieve.getMaxConeSize() < 0) return;
1811:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
1812:       int *orientations = new int[isieve.getMaxConeSize()];

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

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

1821:           orientations[i] = orientation->restrictPoint(arrow)[0];
1822:         }
1823:         isieve.setConeOrientation(orientations, renumbering[*b_iter]);
1824:       }
1825:       delete [] orientations;
1826:     };
1827:     template<typename Section, typename ISection, typename Renumbering>
1828:     static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
1829:       const typename Section::chart_type& chart = coordinates.getChart();
1830:       typename ISection::point_type       min   = *chart.begin();
1831:       typename ISection::point_type       max   = *chart.begin();

1833:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1834:         min = std::min(min, *p_iter);
1835:         max = std::max(max, *p_iter);
1836:       }
1837:       icoordinates.setChart(typename ISection::chart_type(min, max+1));
1838:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1839:         icoordinates.setFiberDimension(*p_iter, coordinates.getFiberDimension(*p_iter));
1840:       }
1841:       icoordinates.allocatePoint();
1842:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1843:         icoordinates.updatePoint(*p_iter, coordinates.restrictPoint(*p_iter));
1844:       }
1845:     };
1846:     template<typename IMesh, typename Label>
1847:     static void convertLabel(IMesh& imesh, const std::string& name, const Obj<Label>& oldLabel) {
1848:       const Obj<typename IMesh::label_type>&        label = imesh.createLabel(name);
1849:       const typename IMesh::sieve_type::chart_type& chart = imesh.getSieve()->getChart();
1850:       int                                           size  = 0;

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

1856:         label->setConeSize(p, coneSize);
1857:         size += coneSize;
1858:       }
1859:       if (size) {label->setSupportSize(0, size);}
1860:       label->allocate();
1861:       for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1862:         const Obj<typename Label::coneSequence>& cone = oldLabel->cone(p);

1864:         if (cone->size()) {
1865:           label->setCone(*cone->begin(), p);
1866:         }
1867:       }
1868:       label->recalculateLabel();
1869:     };
1870:     template<typename Mesh, typename IMesh, typename Renumbering>
1871:     static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
1872:       convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
1873:       imesh.stratify();
1874:       convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
1875:       convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
1876:       const typename Mesh::labels_type& labels = mesh.getLabels();
1877:       std::string heightName("height");
1878:       std::string depthName("depth");

1880:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
1881: #ifdef IMESH_NEW_LABELS
1882:         if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
1883:           convertLabel(imesh, l_iter->first, l_iter->second);
1884:         }
1885: #else
1886:         imesh.setLabel(l_iter->first, l_iter->second);
1887: #endif
1888:       }
1889:     };
1890:   };
1891: }

1893: #endif