Actual source code: ISieve.hh

  1: #ifndef included_ALE_ISieve_hh
  2: #define included_ALE_ISieve_hh

  4: #ifndef  included_ALE_hh
  5: #include <sieve/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:       typedef ISieveVisitor::PointRetriever<Sieve,Retriever> TmpRetriever;
801:       Retriever    pV(200, v, true); // Correct estimate is pow(std::max(1, sieve->getMaxConeSize()), mesh->depth())
802:       TmpRetriever cV[2] = {TmpRetriever(200,pV), TmpRetriever(200,pV)};
803:       int          c     = 0;

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

1298:           for(int c = start; c < end; ++c) {
1299:             this->newCones[q].push_back(this->cones[c]);
1300:           }
1301:         }
1302:         this->newCones[q].push_back(p);
1303:       }
1304:       if (!this->chart.hasPoint(p) || forceSupport) {
1305:         if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1306:           const index_type start = this->supportOffsets[p];
1307:           const index_type end   = this->supportOffsets[p+1];

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

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

1365:         for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1366:           this->cones[c] = oldCones[cO];
1367:         }
1368:         for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1369:           this->supports[s] = oldSupports[sO];
1370:         }
1371:         if (this->orientCones) {
1372:           for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1373:             this->coneOrientations[c] = oldConeOrientations[cO];
1374:           }
1375:         }
1376:       }
1377:       // Inject new cones and supports
1378:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1379:         index_type start = this->coneOffsets[c_iter->first];

1381:         for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1382:           this->cones[start++] = *p_iter;
1383:         }
1384:         if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1385:       }
1386:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1387:         index_type start = this->supportOffsets[s_iter->first];

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

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

1426:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1427:         this->cones[c] = cone[i];
1428:       }
1429:     };
1430:     void setCone(const point_type cone, const point_type& p) {
1431:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1432:       this->chart.checkPoint(p);
1433:       const index_type start = this->coneOffsets[p];
1434:       const index_type end   = this->coneOffsets[p+1];

1436:       if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1437:       this->cones[start] = cone;
1438:     };
1439: #if 0
1440:     template<typename PointSequence>
1441:     void setCone(const PointSequence& cone, const point_type& p) {
1442:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1443:       this->chart.checkPoint(p);
1444:       const index_type start = this->coneOffsets[p];
1445:       const index_type end   = this->coneOffsets[p+1];
1446:       if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1447:       typename PointSequence::iterator c_iter = cone.begin();

1449:       for(index_type c = start; c < end; ++c, ++c_iter) {
1450:         this->cones[c] = *c_iter;
1451:       }
1452:     };
1453: #endif
1454:     void setSupport(const point_type& p, const point_type support[]) {
1455:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1456:       this->chart.checkPoint(p);
1457:       const index_type start = this->supportOffsets[p];
1458:       const index_type end   = this->supportOffsets[p+1];

1460:       for(index_type s = start, i = 0; s < end; ++s, ++i) {
1461:         this->supports[s] = support[i];
1462:       }
1463:     };
1464: #if 0
1465:     template<typename PointSequence>
1466:     void setSupport(const point_type& p, const PointSequence& support) {
1467:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1468:       this->chart.checkPoint(p);
1469:       const index_type start = this->supportOffsets[p];
1470:       const index_type end   = this->supportOffsets[p+1];
1471:       if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1472:       typename PointSequence::iterator s_iter = support.begin();

1474:       for(index_type s = start; s < end; ++s, ++s_iter) {
1475:         this->supports[s] = *s_iter;
1476:       }
1477:     };
1478: #endif
1479:     void setConeOrientation(const int coneOrientation[], const point_type& p) {
1480:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1481:       this->chart.checkPoint(p);
1482:       const index_type start = this->coneOffsets[p];
1483:       const index_type end   = this->coneOffsets[p+1];

1485:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
1486:         this->coneOrientations[c] = coneOrientation[i];
1487:       }
1488:     };
1489:     void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1490:       for(point_type p = 0; p < numCells; ++p) {
1491:         const index_type start = p*numCorners;
1492:         const index_type end   = (p+1)*numCorners;

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

1497:           this->supportOffsets[q+1]++;
1498:         }
1499:       }
1500:       for(point_type p = numCells; p < this->chart.max(); ++p) {
1501:         this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1502:       }
1503:     };
1504:     void symmetrize() {
1505:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1506:       offsets -= this->chart.min();
1507:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1508:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

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

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

1517:           this->chart.checkPoint(q);
1518:           this->supports[this->supportOffsets[q]+offsets[q]] = p;
1519:           ++offsets[q];
1520:         }
1521:       }
1522:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1523:       indexAlloc.deallocate(offsets, this->chart.size()+1);
1524:     }
1525:     index_type getBaseSize() const {
1526:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1527:       return this->baseSize;
1528:     }
1529:     index_type getCapSize() const {
1530:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1531:       return this->capSize;
1532:     }
1533:     template<typename _Section>
1534:     void relabel(_Section& labeling) {

1537:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1538:       offsets -= this->chart.min();
1539:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1540:       // Recalculate coneOffsets
1541:       for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
1542:         const point_type newP = labeling.restrictPoint(p)[0];

1544:         offsets[newP+1] = this->getConeSize(p);
1545:       }
1546:       this->prefixSum(offsets);
1547:       PetscMemcpy(this->coneOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
1548:       // Recalculate supportOffsets
1549:       for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
1550:         const point_type newP = labeling.restrictPoint(p)[0];

1552:         offsets[newP+1] = this->getSupportSize(p);
1553:       }
1554:       this->prefixSum(offsets);
1555:       PetscMemcpy(this->supportOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
1556:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1557:       indexAlloc.deallocate(offsets, this->chart.size()+1);
1558:       index_type  size = std::max(this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()],
1559:                                    this->supportOffsets[this->chart.max()] - this->supportOffsets[this->chart.min()]);
1560:       index_type *orientations = offsets = indexAlloc.allocate(size);
1561:       for(index_type i = 0; i < size; ++i) {indexAlloc.construct(orientations+i, index_type(0));}
1562:       // Recalculate coneOrientations
1563:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1564:         const point_type newP  = labeling.restrictPoint(p)[0];
1565:         const index_type start = this->coneOffsets[newP];
1566:         const index_type end   = this->coneOffsets[newP+1];

1568:         for(index_type c = start; c < end; ++c, ++offset) {
1569:           orientations[c] = this->coneOrientations[offset];
1570:         }
1571:       }
1572:       PetscMemcpy(this->coneOrientations, orientations, (this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()])*sizeof(index_type));CHKERRXX(ierr);
1573:       for(index_type i = 0; i < size; ++i) {indexAlloc.destroy(orientations+i);}
1574:       indexAlloc.deallocate(orientations, size);
1575:       // Recalculate cones
1576:       point_type *array = pointAlloc.allocate(size);

1578:       for(index_type i = 0; i < size; ++i) {pointAlloc.construct(array+i, point_type(0));}
1579:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1580:         const point_type newP  = labeling.restrictPoint(p)[0];
1581:         const index_type start = this->coneOffsets[newP];
1582:         const index_type end   = this->coneOffsets[newP+1];

1584:         for(index_type c = start; c < end; ++c, ++offset) {
1585:           const point_type newQ = labeling.restrictPoint(this->cones[offset])[0];

1587:           array[c] = newQ;
1588:         }
1589:       }
1590:       PetscMemcpy(this->cones, array, size*sizeof(point_type));CHKERRXX(ierr);
1591:       // Recalculate supports
1592:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1593:         const point_type newP  = labeling.restrictPoint(p)[0];
1594:         const index_type start = this->supportOffsets[newP];
1595:         const index_type end   = this->supportOffsets[newP+1];

1597:         for(index_type c = start; c < end; ++c, ++offset) {
1598:           const point_type newQ = labeling.restrictPoint(this->supports[offset])[0];

1600:           array[c] = newQ;
1601:         }
1602:       }
1603:       PetscMemcpy(this->supports, array, size*sizeof(point_type));CHKERRXX(ierr);
1604:       for(index_type i = 0; i < size; ++i) {pointAlloc.destroy(array+i);}
1605:       pointAlloc.deallocate(array, size);
1606:     }
1607:   public: // Traversals
1608:     template<typename Visitor>
1609:     void roots(const Visitor& v) const {
1610:       this->roots(const_cast<Visitor&>(v));
1611:     }
1612:     template<typename Visitor>
1613:     void roots(Visitor& v) const {
1614:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1616:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1617:         if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1618:           if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1619:             v.visitPoint(p);
1620:           }
1621:         }
1622:       }
1623:     }
1624:     template<typename Visitor>
1625:     void leaves(const Visitor& v) const {
1626:       this->leaves(const_cast<Visitor&>(v));
1627:     }
1628:     template<typename Visitor>
1629:     void leaves(Visitor& v) const {
1630:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1632:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1633:         if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1634:           if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1635:             v.visitPoint(p);
1636:           }
1637:         }
1638:       }
1639:     }
1640:     template<typename Visitor>
1641:     void base(const Visitor& v) const {
1642:       this->base(const_cast<Visitor&>(v));
1643:     }
1644:     template<typename Visitor>
1645:     void base(Visitor& v) const {
1646:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1648:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1649:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1650:           v.visitPoint(p);
1651:         }
1652:       }
1653:     }
1654:     template<typename Visitor>
1655:     void cap(const Visitor& v) const {
1656:       this->cap(const_cast<Visitor&>(v));
1657:     }
1658:     template<typename Visitor>
1659:     void cap(Visitor& v) const {
1660:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

1662:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1663:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1664:           v.visitPoint(p);
1665:         }
1666:       }
1667:     }
1668:     template<typename PointSequence, typename Visitor>
1669:     void cone(const PointSequence& points, Visitor& v) const {
1670:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1671:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1672:         const point_type p = *p_iter;
1673:         this->chart.checkPoint(p);
1674:         const index_type start = this->coneOffsets[p];
1675:         const index_type end   = this->coneOffsets[p+1];

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

1690:       for(index_type c = start; c < end; ++c) {
1691:         v.visitArrow(arrow_type(this->cones[c], p));
1692:         v.visitPoint(this->cones[c]);
1693:       }
1694:     }
1695:     const Obj<coneSequence>& cone(const point_type& p) const {
1696:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1697:       if (!this->chart.hasPoint(p)) {
1698:         this->coneSeq->reset(this->cones, 0, 0);
1699:       } else {
1700:         this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1701:       }
1702:       return this->coneSeq;
1703:     }
1704:     template<typename PointSequence, typename Visitor>
1705:     void support(const PointSequence& points, Visitor& v) const {
1706:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1707:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1708:         const point_type p = *p_iter;
1709:         this->chart.checkPoint(p);
1710:         const index_type start = this->supportOffsets[p];
1711:         const index_type end   = this->supportOffsets[p+1];

1713:         for(index_type s = start; s < end; ++s) {
1714:           v.visitArrow(arrow_type(p, this->supports[s]));
1715:           v.visitPoint(this->supports[s]);
1716:         }
1717:       }
1718:     }
1719:     template<typename Visitor>
1720:     void support(const point_type& p, Visitor& v) const {
1721:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1722:       this->chart.checkPoint(p);
1723:       const index_type start = this->supportOffsets[p];
1724:       const index_type end   = this->supportOffsets[p+1];

1726:       for(index_type s = start; s < end; ++s) {
1727:         v.visitArrow(arrow_type(p, this->supports[s]));
1728:         v.visitPoint(this->supports[s]);
1729:       }
1730:     }
1731:     const Obj<supportSequence>& support(const point_type& p) const {
1732:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1733:       if (!this->chart.hasPoint(p)) {
1734:         this->supportSeq->reset(this->supports, 0, 0);
1735:       } else {
1736:         this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1737:       }
1738:       return this->supportSeq;
1739:     }
1740:     template<typename Visitor>
1741:     void orientedCone(const point_type& p, Visitor& v) const {
1742:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1743:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1744:       this->chart.checkPoint(p);
1745:       const index_type start = this->coneOffsets[p];
1746:       const index_type end   = this->coneOffsets[p+1];

1748:       for(index_type c = start; c < end; ++c) {
1749:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1750:         v.visitPoint(this->cones[c], this->coneOrientations[c]);
1751:       }
1752:     }
1753:     template<typename Visitor>
1754:     void orientedReverseCone(const point_type& p, Visitor& v) const {
1755:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1756:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1757:       this->chart.checkPoint(p);
1758:       const index_type start = this->coneOffsets[p];
1759:       const index_type end   = this->coneOffsets[p+1];

1761:       for(index_type c = end-1; c >= start; --c) {
1762:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1763:         v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1764:       }
1765:     }
1766:     template<typename Visitor>
1767:     void orientedSupport(const point_type& p, Visitor& v) const {
1768:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1769:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1770:       this->chart.checkPoint(p);
1771:       const index_type start = this->supportOffsets[p];
1772:       const index_type end   = this->supportOffsets[p+1];

1774:       for(index_type s = start; s < end; ++s) {
1775:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1776:         //v.visitPoint(this->supports[s], this->supportOrientations[s]);
1777:         v.visitArrow(arrow_type(this->supports[s], p), 0);
1778:         v.visitPoint(this->supports[s], 0);
1779:       }
1780:     }
1781:     template<typename Visitor>
1782:     void orientedReverseSupport(const point_type& p, Visitor& v) const {
1783:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1784:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1785:       this->chart.checkPoint(p);
1786:       const index_type start = this->supportOffsets[p];
1787:       const index_type end   = this->supportOffsets[p+1];

1789:       for(index_type s = end-1; s >= start; --s) {
1790:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1791:         //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
1792:         v.visitArrow(arrow_type(this->supports[s], p), 0);
1793:         v.visitPoint(this->supports[s], 0);
1794:       }
1795:     }
1796:     // Currently does only 1 level
1797:     //   Does not check for uniqueness
1798:     template<typename Visitor>
1799:     void meet(const point_type& p, const point_type& q, Visitor& v) const {
1800:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1801:       this->chart.checkPoint(p);
1802:       this->chart.checkPoint(q);
1803:       const index_type startP = this->coneOffsets[p];
1804:       const index_type endP   = this->coneOffsets[p+1];
1805:       const index_type startQ = this->coneOffsets[q];
1806:       const index_type endQ   = this->coneOffsets[q+1];

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

1811:         for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1812:           if (c1 == this->cones[cQ]) v.visitPoint(c1);
1813:         }
1814:         if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1815:       }
1816:     }
1817:     // Currently does only 1 level
1818:     template<typename Sequence, typename Visitor>
1819:     void join(const Sequence& points, Visitor& v) const {
1820:       typedef std::set<point_type> pointSet;
1821:       pointSet intersect[2] = {pointSet(), pointSet()};
1822:       pointSet tmp;
1823:       int      p = 0;
1824:       int      c = 0;

1826:       for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1827:         this->chart.checkPoint(*p_iter);
1828:         tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1829:         if (p == 0) {
1830:           intersect[1-c].insert(tmp.begin(), tmp.end());
1831:           p++;
1832:         } else {
1833:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1834:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1835:           intersect[c].clear();
1836:         }
1837:         c = 1 - c;
1838:         tmp.clear();
1839:       }
1840:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1841:         v.visitPoint(*p_iter);
1842:       }
1843:     }
1844:     // Helper function
1845:     void insertNSupport(point_type p, pointSet& set, const int depth) {
1846:       const index_type start = this->supportOffsets[p];
1847:       const index_type end   = this->supportOffsets[p+1];

1849:       if (depth == 1) {
1850:         set.insert(&this->supports[start], &this->supports[end]);
1851:       } else {
1852:         for(index_type s = start; s < end; ++s) {
1853:           this->insertNSupport(this->supports[s], set, depth-1);
1854:         }
1855:       }
1856:     }
1857:     // Gives only the join of depth n
1858:     template<typename SequenceIterator, typename Visitor>
1859:     void nJoin(const SequenceIterator& pointsBegin, const SequenceIterator& pointsEnd, const int depth, Visitor& v) {
1860:       typedef std::set<point_type> pointSet;
1861:       pointSet intersect[2] = {pointSet(), pointSet()};
1862:       pointSet tmp;
1863:       int      p = 0;
1864:       int      c = 0;

1866:       for(SequenceIterator p_iter = pointsBegin; p_iter != pointsEnd; ++p_iter) {
1867:         this->chart.checkPoint(*p_iter);
1868:         // Put points in the nSupport into tmp (duplicates are fine since it is a set)
1869:         this->insertNSupport(*p_iter, tmp, depth);
1870:         if (p == 0) {
1871:           intersect[1-c].insert(tmp.begin(), tmp.end());
1872:           p++;
1873:         } else {
1874:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1875:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1876:           intersect[c].clear();
1877:         }
1878:         c = 1 - c;
1879:         tmp.clear();
1880:       }
1881:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1882:         v.visitPoint(*p_iter);
1883:       }
1884:     }
1885:   public: // Viewing
1886:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1887:       ostringstream txt;
1888:       int rank;

1890:       if (comm == MPI_COMM_NULL) {
1891:         comm = this->comm();
1892:         rank = this->commRank();
1893:       } else {
1894:         MPI_Comm_rank(comm, &rank);
1895:       }
1896:       if (name == "") {
1897:         if(rank == 0) {
1898:           txt << "viewing an IFSieve" << std::endl;
1899:         }
1900:       } else {
1901:         if(rank == 0) {
1902:           txt << "viewing IFSieve '" << name << "'" << std::endl;
1903:         }
1904:       }
1905:       PetscSynchronizedPrintf(comm, "Max sizes cone: %d support: %d\n", this->getMaxConeSize(), this->getMaxSupportSize());
1906:       if(rank == 0) {
1907:         txt << "cap --> base:" << std::endl;
1908:       }
1909:       ISieveVisitor::PrintVisitor pV(txt, rank);
1910:       this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1911:       PetscSynchronizedPrintf(comm, txt.str().c_str());
1912:       PetscSynchronizedFlush(comm);
1913:       ostringstream txt2;

1915:       if(rank == 0) {
1916:         txt2 << "base <-- cap:" << std::endl;
1917:       }
1918:       ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1919:       this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1920:       PetscSynchronizedPrintf(comm, txt2.str().c_str());
1921:       PetscSynchronizedFlush(comm);
1922:       if (orientCones) {
1923:         ostringstream txt3;

1925:         if(rank == 0) {
1926:           txt3 << "Orientation:" << std::endl;
1927:         }
1928:         ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1929:         this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1930:         PetscSynchronizedPrintf(comm, txt3.str().c_str());
1931:         PetscSynchronizedFlush(comm);
1932:       }
1933:     };
1934:   };

1936:   class ISieveConverter {
1937:   public:
1938:     template<typename Sieve, typename ISieve, typename Renumbering>
1939:     static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1940:       // First construct a renumbering of the sieve points
1941:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
1942:       const Obj<typename Sieve::capSequence>&  cap  = sieve.cap();
1943:       typename ISieve::point_type              min  = 0;
1944:       typename ISieve::point_type              max  = 0;

1946:       if (renumber) {
1947:         /* Roots/Leaves from Sieve do not seem to work */

1949:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1950:           if (sieve.support(*b_iter)->size() == 0) {
1951:             renumbering[*b_iter] = max++;
1952:           }
1953:         }
1954:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1955:           if (sieve.cone(*c_iter)->size() == 0) {
1956:             renumbering[*c_iter] = max++;
1957:           }
1958:         }
1959:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1960:           if (renumbering.find(*b_iter) == renumbering.end()) {
1961:             renumbering[*b_iter] = max++;
1962:           }
1963:         }
1964:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1965:           if (renumbering.find(*c_iter) == renumbering.end()) {
1966:             renumbering[*c_iter] = max++;
1967:           }
1968:         }
1969:       } else {
1970:         if (base->size()) {
1971:           min = *base->begin();
1972:           max = *base->begin();
1973:         } else if (cap->size()) {
1974:           min = *cap->begin();
1975:           max = *cap->begin();
1976:         }
1977:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1978:           min = std::min(min, *b_iter);
1979:           max = std::max(max, *b_iter);
1980:         }
1981:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1982:           min = std::min(min, *c_iter);
1983:           max = std::max(max, *c_iter);
1984:         }
1985:         if (base->size() || cap->size()) {
1986:           ++max;
1987:         }
1988:         for(typename ISieve::point_type p = min; p < max; ++p) {
1989:           renumbering[p] = p;
1990:         }
1991:       }
1992:       // Create the ISieve
1993:       isieve.setChart(typename ISieve::chart_type(min, max));
1994:       // Set cone and support sizes
1995:       size_t maxSize = 0;

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

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

2006:         isieve.setSupportSize(renumbering[*c_iter], support->size());
2007:         maxSize = std::max(maxSize, support->size());
2008:       }
2009:       isieve.allocate();
2010:       // Fill up cones and supports
2011:       typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];

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

2017:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2018:           points[i] = renumbering[*c_iter];
2019:         }
2020:         isieve.setCone(points, renumbering[*b_iter]);
2021:       }
2022:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2023:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2024:         int i = 0;

2026:         for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
2027:           points[i] = renumbering[*s_iter];
2028:         }
2029:         isieve.setSupport(renumbering[*c_iter], points);
2030:       }
2031:       delete [] points;
2032:     }
2033:     template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
2034:     static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
2035:       if (isieve.getMaxConeSize() < 0) return;
2036:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
2037:       int *orientations = new int[isieve.getMaxConeSize()];

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

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

2046:           orientations[i] = orientation->restrictPoint(arrow)[0];
2047:         }
2048:         isieve.setConeOrientation(orientations, renumbering[*b_iter]);
2049:       }
2050:       delete [] orientations;
2051:     }
2052:     template<typename Section, typename ISection, typename Renumbering>
2053:     static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
2054:       const typename Section::chart_type& chart = coordinates.getChart();
2055:       typename ISection::point_type       min   = *chart.begin();
2056:       typename ISection::point_type       max   = *chart.begin();

2058:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2059:         min = std::min(min, renumbering[*p_iter]);
2060:         max = std::max(max, renumbering[*p_iter]);
2061:       }
2062:       icoordinates.setChart(typename ISection::chart_type(min, max+1));
2063:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2064:         icoordinates.setFiberDimension(renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));
2065:       }
2066:       icoordinates.allocatePoint();
2067:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2068:         icoordinates.updatePoint(renumbering[*p_iter], coordinates.restrictPoint(*p_iter));
2069:       }
2070:     }
2071:     template<typename Label, typename Renumbering>
2072:     static void convertLabel(const Obj<Label>& newLabel, const Obj<Label>& oldLabel, Renumbering& renumbering) {
2073:       for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
2074:         if (oldLabel->getConeSize(p->first)) {
2075:           newLabel->setCone(*oldLabel->cone(p->first)->begin(), p->second);
2076:         }
2077:       }
2078:     }
2079:     template<typename Mesh, typename IMesh, typename Renumbering>
2080:     static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
2081:       convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
2082:       imesh.stratify();
2083:       convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
2084:       convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
2085:       if (mesh.hasRealSection("normals")) {
2086:         convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
2087:       }
2088:       const typename Mesh::labels_type& labels = mesh.getLabels();
2089:       std::string heightName("height");
2090:       std::string depthName("depth");

2092:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
2093: #ifdef IMESH_NEW_LABELS
2094:         if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
2095:           convertLabel(imesh, l_iter->first, l_iter->second);
2096:         }
2097: #else
2098:         if (renumber) {
2099:           convertLabel(imesh.createLabel(l_iter->first), l_iter->second, renumbering);
2100:         } else {
2101:           imesh.setLabel(l_iter->first, l_iter->second);
2102:         }
2103: #endif
2104:       }
2105:     }
2106:   };

2108:   class ISieveSerializer {
2109:   public:
2110:     template<typename ISieve>
2111:     static void writeSieve(const std::string& filename, ISieve& sieve) {
2112:       std::ofstream fs;

2114:       if (sieve.commRank() == 0) {
2115:         fs.open(filename.c_str());
2116:       }
2117:       writeSieve(fs, sieve);
2118:       if (sieve.commRank() == 0) {
2119:         fs.close();
2120:       }
2121:     }
2122:     template<typename ISieve>
2123:     static void writeSieve(std::ofstream& fs, ISieve& sieve) {
2124:       typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
2125:       const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
2126:       typename ISieve::point_type             min   = chart->min();
2127:       typename ISieve::point_type             max   = chart->max();
2128:       PetscInt                               *mins  = new PetscInt[sieve.commSize()];
2129:       PetscInt                               *maxs  = new PetscInt[sieve.commSize()];

2131:       // Write sizes
2132:       if (sieve.commRank() == 0) {
2133:         // Write local
2134:         fs << min <<" "<< max << std::endl;
2135:         for(typename ISieve::point_type p = min; p < max; ++p) {
2136:           fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
2137:         }
2138:         // Receive and write remote
2139:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2140:           PetscInt       minp, maxp;
2141:           PetscInt      *sizes;
2142:           MPI_Status     status;

2145:           MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2146:           MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2147:           PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2148:           MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2149:           fs << minp <<" "<< maxp << std::endl;
2150:           for(PetscInt s = 0; s < maxp - minp; ++s) {
2151:             fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
2152:           }
2153:           PetscFree(sizes);CHKERRXX(ierr);
2154:           mins[pr] = minp;
2155:           maxs[pr] = maxp;
2156:         }
2157:       } else {
2158:         // Send remote
2159:         //PetscInt       min = chart->min();
2160:         //PetscInt       max = chart->max();
2161:         PetscInt       s   = 0;
2162:         PetscInt      *sizes;

2165:         MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2166:         MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2167:         PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
2168:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2169:           sizes[s*2+0] = sieve.getConeSize(p);
2170:           sizes[s*2+1] = sieve.getSupportSize(p);
2171:         }
2172:         MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2173:         PetscFree(sizes);CHKERRXX(ierr);
2174:       }
2175:       // Write covering
2176:       if (sieve.commRank() == 0) {
2177:         // Write local
2178:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);

2180:         for(typename ISieve::point_type p = min; p < max; ++p) {
2181:           sieve.cone(p, pV);
2182:           const typename Visitor::point_type *cone  = pV.getPoints();
2183:           const int                           cSize = pV.getSize();

2185:           if (cSize > 0) {
2186:             for(int c = 0; c < cSize; ++c) {
2187:               if (c) {fs << " ";}
2188:               fs << cone[c];
2189:             }
2190:             fs << std::endl;
2191:           }
2192:           pV.clear();

2194:           sieve.orientedCone(p, pV);
2195:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2196:           const int                                    oSize = pV.getOrientedSize();

2198:           if (oSize > 0) {
2199:             for(int c = 0; c < oSize; ++c) {
2200:               if (c) {fs << " ";}
2201:               fs << oCone[c].second;
2202:             }
2203:             fs << std::endl;
2204:           }
2205:           pV.clear();

2207:           sieve.support(p, pV);
2208:           const typename Visitor::point_type *support = pV.getPoints();
2209:           const int                           sSize   = pV.getSize();

2211:           if (sSize > 0) {
2212:             for(int s = 0; s < sSize; ++s) {
2213:               if (s) {fs << " ";}
2214:               fs << support[s];
2215:             }
2216:             fs << std::endl;
2217:           }
2218:           pV.clear();
2219:         }
2220:         // Receive and write remote
2221:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2222:           PetscInt       size;
2223:           PetscInt      *data;
2224:           PetscInt       off = 0;
2225:           MPI_Status     status;

2228:           MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2229:           PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2230:           MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2231:           for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2232:             PetscInt cSize = data[off++];

2234:             fs << cSize << std::endl;
2235:             if (cSize > 0) {
2236:               for(int c = 0; c < cSize; ++c) {
2237:                 if (c) {fs << " ";}
2238:                 fs << data[off++];
2239:               }
2240:               fs << std::endl;
2241:             }
2242:             PetscInt oSize = data[off++];

2244:             if (oSize > 0) {
2245:               for(int c = 0; c < oSize; ++c) {
2246:                 if (c) {fs << " ";}
2247:                 fs << data[off++];
2248:               }
2249:               fs << std::endl;
2250:             }
2251:             PetscInt sSize = data[off++];

2253:             fs << sSize << std::endl;
2254:             if (sSize > 0) {
2255:               for(int s = 0; s < sSize; ++s) {
2256:                 if (s) {fs << " ";}
2257:                 fs << data[off++];
2258:               }
2259:               fs << std::endl;
2260:             }
2261:           }
2262:           assert(off == size);
2263:           PetscFree(data);CHKERRXX(ierr);
2264:         }
2265:       } else {
2266:         // Send remote
2267:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2268:         PetscInt totalConeSize    = 0;
2269:         PetscInt totalSupportSize = 0;

2271:         for(typename ISieve::point_type p = min; p < max; ++p) {
2272:           totalConeSize    += sieve.getConeSize(p);
2273:           totalSupportSize += sieve.getSupportSize(p);
2274:         }
2275:         PetscInt       size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
2276:         PetscInt       off  = 0;
2277:         PetscInt      *data;

2280:         MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2281:         // There is no nice way to make a generic MPI type here. Really sucky
2282:         PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2283:         for(typename ISieve::point_type p = min; p < max; ++p) {
2284:           sieve.cone(p, pV);
2285:           const typename Visitor::point_type *cone  = pV.getPoints();
2286:           const int                           cSize = pV.getSize();

2288:           data[off++] = cSize;
2289:           for(int c = 0; c < cSize; ++c) {
2290:             data[off++] = cone[c];
2291:           }
2292:           pV.clear();

2294:           sieve.orientedCone(p, pV);
2295:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2296:           const int                                    oSize = pV.getOrientedSize();

2298:           data[off++] = oSize;
2299:           for(int c = 0; c < oSize; ++c) {
2300:             data[off++] = oCone[c].second;
2301:           }
2302:           pV.clear();

2304:           sieve.support(p, pV);
2305:           const typename Visitor::point_type *support = pV.getPoints();
2306:           const int                           sSize   = pV.getSize();

2308:           data[off++] = sSize;
2309:           for(int s = 0; s < sSize; ++s) {
2310:             data[off++] = support[s];
2311:           }
2312:           pV.clear();
2313:         }
2314:         MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2315:         PetscFree(data);CHKERRXX(ierr);
2316:       }
2317:       delete [] mins;
2318:       delete [] maxs;
2319:       // Output renumbering
2320:     }
2321:     template<typename ISieve>
2322:     static void loadSieve(const std::string& filename, ISieve& sieve) {
2323:       std::ifstream fs;

2325:       if (sieve.commRank() == 0) {
2326:         fs.open(filename.c_str());
2327:       }
2328:       loadSieve(fs, sieve);
2329:       if (sieve.commRank() == 0) {
2330:         fs.close();
2331:       }
2332:     }
2333:     template<typename ISieve>
2334:     static void loadSieve(std::ifstream& fs, ISieve& sieve) {
2335:       typename ISieve::point_type min, max;
2336:       PetscInt                   *mins = new PetscInt[sieve.commSize()];
2337:       PetscInt                   *maxs = new PetscInt[sieve.commSize()];
2338:       PetscInt                   *totalConeSizes    = new PetscInt[sieve.commSize()];
2339:       PetscInt                   *totalSupportSizes = new PetscInt[sieve.commSize()];

2341:       // Load sizes
2342:       if (sieve.commRank() == 0) {
2343:         // Load local
2344:         fs >> min;
2345:         fs >> max;
2346:         sieve.setChart(typename ISieve::chart_type(min, max));
2347:         for(typename ISieve::point_type p = min; p < max; ++p) {
2348:           typename ISieve::index_type coneSize, supportSize;

2350:           fs >> coneSize;
2351:           fs >> supportSize;
2352:           sieve.setConeSize(p, coneSize);
2353:           sieve.setSupportSize(p, supportSize);
2354:         }
2355:         // Load and send remote
2356:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2357:           PetscInt       minp, maxp;
2358:           PetscInt      *sizes;

2361:           fs >> minp;
2362:           fs >> maxp;
2363:           MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2364:           MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2365:           PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2366:           totalConeSizes[pr]    = 0;
2367:           totalSupportSizes[pr] = 0;
2368:           for(PetscInt s = 0; s < maxp - minp; ++s) {
2369:             fs >> sizes[s*2+0];
2370:             fs >> sizes[s*2+1];
2371:             totalConeSizes[pr]    += sizes[s*2+0];
2372:             totalSupportSizes[pr] += sizes[s*2+1];
2373:           }
2374:           MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2375:           PetscFree(sizes);CHKERRXX(ierr);
2376:           mins[pr] = minp;
2377:           maxs[pr] = maxp;
2378:         }
2379:       } else {
2380:         // Load remote
2381:         PetscInt       s   = 0;
2382:         PetscInt      *sizes;
2383:         MPI_Status     status;

2386:         MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2387:         MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2388:         sieve.setChart(typename ISieve::chart_type(min, max));
2389:         PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2390:         MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2391:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2392:           sieve.setConeSize(p, sizes[s*2+0]);
2393:           sieve.setSupportSize(p, sizes[s*2+1]);
2394:         }
2395:         PetscFree(sizes);CHKERRXX(ierr);
2396:       }
2397:       sieve.allocate();
2398:       // Load covering
2399:       if (sieve.commRank() == 0) {
2400:         // Load local
2401:         typename ISieve::index_type  maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
2402:         typename ISieve::point_type *points  = new typename ISieve::point_type[maxSize];

2404:         for(typename ISieve::point_type p = min; p < max; ++p) {
2405:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
2406:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);

2408:           if (coneSize > 0) {
2409:             for(int c = 0; c < coneSize; ++c) {
2410:               fs >> points[c];
2411:             }
2412:             sieve.setCone(points, p);
2413:             if (sieve.orientedCones()) {
2414:               for(int c = 0; c < coneSize; ++c) {
2415:                 fs >> points[c];
2416:               }
2417:               sieve.setConeOrientation(points, p);
2418:             }
2419:           }
2420:           if (supportSize > 0) {
2421:             for(int s = 0; s < supportSize; ++s) {
2422:               fs >> points[s];
2423:             }
2424:             sieve.setSupport(p, points);
2425:           }
2426:         }
2427:         delete [] points;
2428:         // Load and send remote
2429:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
2430:           PetscInt       size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
2431:           PetscInt       off  = 0;
2432:           PetscInt      *data;

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

2441:             fs >> coneSize;
2442:             data[off++] = coneSize;
2443:             if (coneSize > 0) {
2444:               for(int c = 0; c < coneSize; ++c) {
2445:                 fs >> data[off++];
2446:               }
2447:               for(int c = 0; c < coneSize; ++c) {
2448:                 fs >> data[off++];
2449:               }
2450:             }
2451:             fs >> supportSize;
2452:             data[off++] = supportSize;
2453:             if (supportSize > 0) {
2454:               for(int s = 0; s < supportSize; ++s) {
2455:                 fs >> data[off++];
2456:               }
2457:             }
2458:           }
2459:           assert(off == size);
2460:           MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2461:           PetscFree(data);CHKERRXX(ierr);
2462:         }
2463:         delete [] mins;
2464:         delete [] maxs;
2465:       } else {
2466:         // Load remote
2467:         PetscInt       size;
2468:         PetscInt      *data;
2469:         PetscInt       off = 0;
2470:         MPI_Status     status;

2473:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2474:         PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2475:         MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2476:         for(typename ISieve::point_type p = min; p < max; ++p) {
2477:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
2478:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2479:           PetscInt cs = data[off++];

2481:           assert(cs == coneSize);
2482:           if (coneSize > 0) {
2483:             sieve.setCone(&data[off], p);
2484:             off += coneSize;
2485:             if (sieve.orientedCones()) {
2486:               sieve.setConeOrientation(&data[off], p);
2487:               off += coneSize;
2488:             }
2489:           }
2490:           PetscInt ss = data[off++];

2492:           assert(ss == supportSize);
2493:           if (supportSize > 0) {
2494:             sieve.setSupport(p, &data[off]);
2495:             off += supportSize;
2496:           }
2497:         }
2498:         assert(off == size);
2499:       }
2500:       // Load renumbering
2501:     }
2502:   };
2503: }

2505: #endif