Actual source code: ISieve.hh
petsc-3.4.2 2013-07-02

      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 <petscdmplex.h>

     10: #include <fstream>

     12: //#define IMESH_NEW_LABELS

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

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

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

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

    570:         if (orientation >= 0) {
    571:           for(int d = 0; d < dim; ++d, ++i) {
    572:             this->values[i] = v[d];
    573:           }
    574:         } else {
    575:           for(int d = dim-1; d >= 0; --d, ++i) {
    576:             this->values[i] = v[d];
    577:           }
    578:         }
    579:       }
    580:       template<typename Arrow>
    581:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    582:     public:
    583:       const value_type *getValues() const {return this->values;};
    584:       int  getSize() const {return this->i;};
    585:       int  getMaxSize() const {return this->size;};
    586:       void ensureSize(const int size) {
    587:         this->clear();
    588:         if (size > this->size) {
    589:           this->size = size;
    590:           if (this->allocated) {delete [] this->values;}
    591:           this->values = new value_type[this->size];
    592:           this->allocated = true;
    593:         }
    594:       };
    595:       void clear() {this->i = 0;};
    596:     };
    597:     template<typename ValueType>
    598:     class RestrictVecVisitor {
    599:     public:
    600:       typedef ValueType value_type;
    601:     protected:
    602:       const Vec          v;
    603:       const PetscSection section;
    604:       int                size;
    605:       int                i;
    606:       int                nF;
    607:       int               *offsets;
    608:       int               *indices;
    609:       bool               processed;
    610:       value_type        *values;
    611:       bool               allocated;
    612:       value_type        *array;
    613:     protected:
    614:       inline void swap(value_type& v, value_type& w) {
    615:         value_type tmp = v;
    616:         v = w;
    617:         w = tmp;
    618:       };
    619:       void processArray() {
    620:         for(PetscInt f = 1; f < nF; ++f) {
    621:           offsets[f+1] += offsets[f];
    622:         }
    623:         for(PetscInt i = 0; i < offsets[nF]; ++i) {
    624:           indices[i] = offsets[indices[i]]++;
    625:         }
    626:         assert(offsets[nF-1] == offsets[nF]);
    627:         for(PetscInt i = 0; i < offsets[nF]; ++i) {
    628:           if (indices[i] == -1) continue;
    629:           PetscInt   startPos = indices[i];
    630:           PetscInt   j        = startPos, k;
    631:           value_type val      = values[i];

    633:           do {
    634:             swap(val, values[j]);
    635:             k = indices[j];
    636:             indices[j] = -1;
    637:             j = k;
    638:           } while(j != startPos);
    639:         }
    640:       };
    641:     public:
    642:       RestrictVecVisitor(const Vec v, const PetscSection s, const int size) : v(v), section(s), size(size), i(0), nF(0), processed(true) {
    643:         this->values    = new value_type[this->size];
    644:         this->allocated = true;
    645:         PetscErrorCode VecGetArray(this->v, &this->array);CHKERRXX(ierr);
    646:       };
    647:       RestrictVecVisitor(const Vec v, const PetscSection s, const int size, value_type *values) : v(v), section(s), size(size), i(0), nF(0), processed(true) {
    648:         this->values    = values;
    649:         this->allocated = false;
    650:         PetscErrorCode VecGetArray(this->v, &this->array);CHKERRXX(ierr);
    651:       };
    652:       RestrictVecVisitor(const Vec v, const PetscSection s, const int size, value_type *values, int *offsets, int *indices) : v(v), section(s), size(size), i(0), nF(0), offsets(offsets), indices(indices), processed(false) {

    655:         this->values    = values;
    656:         this->allocated = false;
    657:         VecGetArray(this->v, &this->array);CHKERRXX(ierr);
    658:         PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
    659:         for(PetscInt f = 0; f <= nF; ++f) {offsets[f] = 0;}
    660:       };
    661:       ~RestrictVecVisitor() {
    662:         if (this->allocated) {delete [] this->values;}
    663:         PetscErrorCode VecRestoreArray(this->v, &this->array);CHKERRXX(ierr);
    664:       };
    665:       template<typename Point>
    666:       inline void visitPoint(const Point& point, const int orientation) {
    667:         // Known:
    668:         //   Nf: number of fields
    669:         // Unknown:
    670:         //   Np: number of points
    671:         //   P:  points
    672:         // Algorithm:
    673:         //   Pass 1: Stack up values as before, but also
    674:         //           count size of each field
    675:         //   Comp 1: Sum up sizes to get field offsets
    676:         //   Pass 2: Number each entry with its intended position
    677:         //   Pass 3: Reorder entries
    678:         // Algorithm if field sizes are known:
    679:         //   Comp 1: Partition array into field components
    680:         //   Pass 1: Stack up values at field offsets
    681:         PetscInt       dim, off;

    684:         PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
    685:         if (i+dim > size) {
    686:           ostringstream msg;
    687:           msg << "Too many values for RestrictVisitor "<<i+dim<<" > "<<size<< std::endl;
    688:           throw ALE::Exception(msg.str().c_str());
    689:         }
    690:         PetscSectionGetOffset(section, point, &off);CHKERRXX(ierr);
    691:         const value_type *v = &array[off];

    693:         if (nF) {
    694:           for(PetscInt f = 0, fOff = 0; f < nF; ++f) {
    695:             PetscInt comp, fDim;

    697:             PetscSectionGetFieldDof(section, point, f, &fDim);CHKERRXX(ierr);
    698:             offsets[f+1] += fDim;
    699:             for(PetscInt d = 0; d < fDim; ++d) {
    700:               indices[i+d] = f;
    701:             }
    702:             if (orientation >= 0) {
    703:               for(PetscInt d = 0; d < fDim; ++d, ++i) {
    704:                 this->values[i] = v[fOff+d];
    705:               }
    706:             } else {
    707:               PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
    708:               for(PetscInt d = fDim/comp-1; d >= 0; --d) {
    709:                 for(PetscInt c = 0; c < comp; ++c, ++i) {
    710:                   this->values[i] = v[fOff+d*comp+c];
    711:                 }
    712:               }
    713:             }
    714:             fOff += fDim;
    715:           }
    716:         } else {
    717:           if (orientation >= 0) {
    718:             for(PetscInt d = 0; d < dim; ++d, ++i) {
    719:               this->values[i] = v[d];
    720:             }
    721:           } else {
    722:             for(PetscInt d = dim-1; d >= 0; --d, ++i) {
    723:               this->values[i] = v[d];
    724:             }
    725:           }
    726:         }
    727:       }
    728:       template<typename Arrow>
    729:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    730:     public:
    731:       const value_type *getValues() {
    732:         if (!processed) {processArray(); processed = true;}
    733:         return this->values;
    734:       };
    735:       int  getSize() const {return this->i;};
    736:       int  getMaxSize() const {return this->size;};
    737:       void ensureSize(const int size) {
    738:         this->clear();
    739:         if (size > this->size) {
    740:           this->size = size;
    741:           if (this->allocated) {delete [] this->values;}
    742:           this->values = new value_type[this->size];
    743:           this->allocated = true;
    744:         }
    745:       };
    746:       void clear() {this->i = 0; if (processed) {processed = false;}};
    747:     };
    748:     template<typename Section>
    749:     class UpdateVisitor {
    750:     public:
    751:       typedef typename Section::value_type value_type;
    752:     protected:
    753:       Section&          section;
    754:       const value_type *values;
    755:       int               i;
    756:     public:
    757:       UpdateVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
    758:       template<typename Point>
    759:       inline void visitPoint(const Point& point, const int orientation) {
    760:         const int dim = section.getFiberDimension(point);
    761:         this->section.updatePoint(point, &this->values[this->i], orientation);
    762:         this->i += dim;
    763:       }
    764:       template<typename Arrow>
    765:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    766:       void clear() {this->i = 0;};
    767:     };
    768:     template<typename Section>
    769:     class UpdateAllVisitor {
    770:     public:
    771:       typedef typename Section::value_type value_type;
    772:     protected:
    773:       Section&          section;
    774:       const value_type *values;
    775:       int               i;
    776:     public:
    777:       UpdateAllVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
    778:       template<typename Point>
    779:       inline void visitPoint(const Point& point, const int orientation) {
    780:         const int dim = section.getFiberDimension(point);
    781:         this->section.updatePointAll(point, &this->values[this->i], orientation);
    782:         this->i += dim;
    783:       }
    784:       template<typename Arrow>
    785:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    786:       void clear() {this->i = 0;};
    787:     };
    788:     template<typename Section>
    789:     class UpdateAddVisitor {
    790:     public:
    791:       typedef typename Section::value_type value_type;
    792:     protected:
    793:       Section&          section;
    794:       const value_type *values;
    795:       int               i;
    796:     public:
    797:       UpdateAddVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
    798:       template<typename Point>
    799:       inline void visitPoint(const Point& point, const int orientation) {
    800:         const int dim = section.getFiberDimension(point);
    801:         this->section.updateAddPoint(point, &this->values[this->i], orientation);
    802:         this->i += dim;
    803:       }
    804:       template<typename Arrow>
    805:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    806:       void clear() {this->i = 0;};
    807:     };
    808:     template<typename ValueType>
    809:     class UpdateVecVisitor {
    810:     public:
    811:       typedef ValueType value_type;
    812:     protected:
    813:       const Vec          v;
    814:       const PetscSection section;
    815:       const value_type  *values;
    816:       const InsertMode   mode;
    817:       PetscInt           nF;
    818:       PetscInt           i;
    819:       value_type        *array;
    820:       PetscInt          *fieldSize;
    821:       PetscInt          *j;
    822:     protected:
    823:       inline static void add   (value_type& x, value_type y) {x += y;}
    824:       inline static void insert(value_type& x, value_type y) {x  = y;}
    825:       template<typename Point>
    826:       void updatePoint(const Point& point, void (*fuse)(value_type&, value_type), const bool setBC, const int orientation = 1) {
    827:         PetscInt        dim;  // The number of dof on this point
    828:         PetscInt        cDim; // The nubmer of constraints on this point
    829:         const PetscInt *cDof; // The indices of the constrained dofs on this point
    830:         value_type     *a;    // The values on this point
    831:         PetscInt        offset, cInd = 0;
    832:         PetscErrorCode  ierr;

    834:         PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
    835:         PetscSectionGetConstraintDof(section, point, &cDim);CHKERRXX(ierr);
    836:         PetscSectionGetOffset(section, point, &offset);CHKERRXX(ierr);
    837:         a    = &array[offset];
    838:         if (!cDim || setBC) {
    839:           if (orientation >= 0) {
    840:             for(PetscInt k = 0; k < dim; ++k) {
    841:               fuse(a[k], values[i+k]);
    842:             }
    843:           } else {
    844:             for(PetscInt k = 0; k < dim; ++k) {
    845:               fuse(a[k], values[i+dim-k-1]);
    846:             }
    847:           }
    848:         } else {
    849:           PetscSectionGetConstraintIndices(section, point, &cDof);CHKERRXX(ierr);
    850:           if (orientation >= 0) {
    851:             for(PetscInt k = 0; k < dim; ++k) {
    852:               if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
    853:               fuse(a[k], values[i+k]);
    854:             }
    855:           } else {
    856:             for(PetscInt k = 0; k < dim; ++k) {
    857:               if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
    858:               fuse(a[k], values[i+dim-k-1]);
    859:             }
    860:           }
    861:         }
    862:         i += dim;
    863:       }
    864:       template<typename Point>
    865:       void updatePointFields(const Point& point, void (*fuse)(value_type&, value_type), const bool setBC, const int orientation = 1) {
    866:         value_type    *a;
    867:         PetscInt       offset;
    868:         PetscInt       fOff = 0;

    871:         PetscSectionGetOffset(section, point, &offset);CHKERRXX(ierr);
    872:         a    = &array[offset];
    873:         for(PetscInt f = 0; f < nF; ++f) {
    874:           PetscInt    dim;  // The number of dof for field f on this point
    875:           PetscInt    comp; // The number of components for field f on this point
    876:           PetscInt    cDim; // The nubmer of constraints for field f on this point
    877:           const PetscInt *cDof; // The indices of the constrained dofs for field f on this point
    878:           PetscInt    cInd = 0;

    880:           PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
    881:           PetscSectionGetFieldDof(section, point, f, &dim);CHKERRXX(ierr);
    882:           PetscSectionGetFieldConstraintDof(section, point, f, &cDim);CHKERRXX(ierr);
    883:           if (!cDim || setBC) {
    884:             if (orientation >= 0) {
    885:               for(PetscInt k = 0; k < dim; ++k) {
    886:                 fuse(a[fOff+k], values[j[f]+k]);
    887:               }
    888:             } else {
    889:               for(PetscInt k = dim/comp-1; k >= 0; --k) {
    890:                 for(PetscInt c = 0; c < comp; ++c) {
    891:                   fuse(a[fOff+(dim/comp-1-k)*comp+c], values[j[f]+k*comp+c]);
    892:                 }
    893:               }
    894:             }
    895:           } else {
    896:             PetscSectionGetFieldConstraintIndices(section, point, f, &cDof);CHKERRXX(ierr);
    897:             if (orientation >= 0) {
    898:               for(PetscInt k = 0; k < dim; ++k) {
    899:                 if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
    900:                 fuse(a[fOff+k], values[j[f]+k]);
    901:               }
    902:             } else {
    903:               for(PetscInt k = dim/comp-1; k >= 0; --k) {
    904:                 for(PetscInt c = 0; c < comp; ++c) {
    905:                   PetscInt ind = k*comp+c;
    906:                   if ((cInd < cDim) && (ind == cDof[cInd])) {++cInd; continue;}
    907:                   fuse(a[fOff+(dim/comp-1-k)*comp+c], values[j[f]+ind]);
    908:                 }
    909:               }
    910:             }
    911:           }
    912:           fOff += dim;
    913:           j[f] += dim;
    914:         }
    915:       }
    916:     public:
    917:       UpdateVecVisitor(const Vec v, const PetscSection s, const value_type *values, InsertMode mode) : v(v), section(s), values(values), mode(mode), nF(0), i(0) {};
    918:       UpdateVecVisitor(const Vec v, const PetscSection s, const value_type *values, InsertMode mode, PetscInt numFields, PetscInt fieldSize[]) : v(v), section(s), values(values), mode(mode), nF(numFields), i(0) {

    921:         VecGetArray(this->v, &this->array);CHKERRXX(ierr);
    922:         PetscMalloc2(numFields,PetscInt,&this->fieldSize,numFields,PetscInt,&j);CHKERRXX(ierr);
    923:         for(PetscInt f = 0; f < nF; ++f) {
    924:           this->fieldSize[f] = fieldSize[f];
    925:         }
    926:         this->clear();
    927:       };
    928:       ~UpdateVecVisitor() {
    930:         VecRestoreArray(this->v, &this->array);CHKERRXX(ierr);
    931:         PetscFree2(fieldSize,j);CHKERRXX(ierr);
    932:       };
    933:       template<typename Point>
    934:       inline void visitPoint(const Point& point, const int orientation) {
    935:         if (nF) {
    936:           switch(mode) {
    937:           case INSERT_VALUES:
    938:             updatePointFields(point, this->insert, false, orientation);break;
    939:           case INSERT_ALL_VALUES:
    940:             updatePointFields(point, this->insert, true,  orientation);break;
    941:           case ADD_VALUES:
    942:             updatePointFields(point, this->add, false, orientation);break;
    943:           case ADD_ALL_VALUES:
    944:             updatePointFields(point, this->add, true,  orientation);break;
    945:           default:
    946:             throw PETSc::Exception("Invalid mode");
    947:           }
    948:         } else {
    949:           switch(mode) {
    950:           case INSERT_VALUES:
    951:             updatePoint(point, this->insert, false, orientation);break;
    952:           case INSERT_ALL_VALUES:
    953:             updatePoint(point, this->insert, true,  orientation);break;
    954:           case ADD_VALUES:
    955:             updatePoint(point, this->add, false, orientation);break;
    956:           case ADD_ALL_VALUES:
    957:             updatePoint(point, this->add, true,  orientation);break;
    958:           default:
    959:             throw PETSc::Exception("Invalid mode");
    960:           }
    961:         }
    962:       }
    963:       template<typename Arrow>
    964:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    965:     public:
    966:       void clear() {
    967:         this->i = 0;
    968:         if (nF) {
    969:           j[0] = 0;
    970:           for(PetscInt f = 1; f < nF; ++f) {
    971:             j[f] = j[f-1] + fieldSize[f-1];
    972:           }
    973:         }
    974:       };
    975:     };
    976:     template<typename Section, typename Order, typename Value>
    977:     class IndicesVisitor {
    978:     public:
    979:       typedef Value                        value_type;
    980:       typedef typename Section::point_type point_type;
    981:     protected:
    982:       const Section& section;
    983:       // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
    984:       Order&         order;
    985:       int            size;
    986:       int            i, p;
    987:       // If false, constrained indices are returned as negative values. Otherwise, they are omitted
    988:       bool           freeOnly;
    989:       // If true, it allows space for constrained variables (even if the indices are not returned) Wierd
    990:       bool           skipConstraints;
    991:       value_type    *values;
    992:       bool           allocated;
    993:       point_type    *points;
    994:       bool           allocatedPoints;
    995:     protected:
    996:       void getUnconstrainedIndices(const point_type& p, const int orientation) {
    997:         if (i+section.getFiberDimension(p) > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
    998:         if (orientation >= 0) {
    999:           const int start = this->order.getIndex(p);
    1000:           const int end   = start + section.getFiberDimension(p);

    1002:           for(int j = start; j < end; ++j, ++i) {
    1003:             this->values[i] = j;
    1004:           }
    1005:         } else if (!section.getNumSpaces()) {
    1006:           const int start = this->order.getIndex(p);
    1007:           const int end   = start + section.getFiberDimension(p);

    1009:           for(int j = end-1; j >= start; --j, ++i) {
    1010:             this->values[i] = j;
    1011:           }
    1012:         } else {
    1013:           const int numSpaces = section.getNumSpaces();
    1014:           int       start     = this->order.getIndex(p);

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

    1019:             for(int j = end-1; j >= start; --j, ++i) {
    1020:               this->values[i] = j;
    1021:             }
    1022:             start = end;
    1023:           }
    1024:         }
    1025:       };
    1026:       void getConstrainedIndices(const point_type& p, const int orientation) {
    1027:         const int cDim = this->section.getConstraintDimension(p);
    1028:         if (i+cDim > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
    1029:         typedef typename Section::bc_type::value_type index_type;
    1030:         const index_type *cDof  = this->section.getConstraintDof(p);
    1031:         const int         start = this->order.getIndex(p);

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

    1036:           for(int j = start, cInd = 0, k = 0; k < dim; ++k) {
    1037:             if ((cInd < cDim) && (k == cDof[cInd])) {
    1038:               if (!freeOnly) values[i++] = -(k+1);
    1039:               if (skipConstraints) ++j;
    1040:               ++cInd;
    1041:             } else {
    1042:               values[i++] = j++;
    1043:             }
    1044:           }
    1045:         } else {
    1046:           int offset  = 0;
    1047:           int cOffset = 0;
    1048:           int k       = -1;

    1050:           for(int space = 0; space < section.getNumSpaces(); ++space) {
    1051:             const int  dim = this->section.getFiberDimension(p, space);
    1052:             const int tDim = this->section.getConstrainedFiberDimension(p, space);
    1053:             int       cInd = (dim - tDim)-1;

    1055:             k += dim;
    1056:             for(int l = 0, j = start+tDim+offset; l < dim; ++l, --k) {
    1057:               if ((cInd >= 0) && (k == cDof[cInd+cOffset])) {
    1058:                 if (!freeOnly) values[i++] = -(offset+l+1);
    1059:                 if (skipConstraints) --j;
    1060:                 --cInd;
    1061:               } else {
    1062:                 values[i++] = --j;
    1063:               }
    1064:             }
    1065:             k       += dim;
    1066:             offset  += dim;
    1067:             cOffset += dim - tDim;
    1068:           }
    1069:         }
    1070:       };
    1071:     public:
    1072:       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) {
    1073:         this->values    = new value_type[this->size];
    1074:         this->allocated = true;
    1075:         if (unique) {
    1076:           this->points          = new point_type[this->size];
    1077:           this->allocatedPoints = true;
    1078:         } else {
    1079:           this->points          = NULL;
    1080:           this->allocatedPoints = false;
    1081:         }
    1082:       };
    1083:       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) {
    1084:         this->values    = values;
    1085:         this->allocated = false;
    1086:         if (unique) {
    1087:           this->points          = new point_type[this->size];
    1088:           this->allocatedPoints = true;
    1089:         } else {
    1090:           this->points          = NULL;
    1091:           this->allocatedPoints = false;
    1092:         }
    1093:       };
    1094:       ~IndicesVisitor() {
    1095:         if (this->allocated) {delete [] this->values;}
    1096:         if (this->allocatedPoints) {delete [] this->points;}
    1097:       };
    1098:       inline void visitPoint(const point_type& point, const int orientation) {
    1099:         if (p >= size) {
    1100:           ostringstream msg;
    1101:           msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
    1102:           throw ALE::Exception(msg.str().c_str());
    1103:         }
    1104:         if (points) {
    1105:           int pp;
    1106:           for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
    1107:           if (pp != p) return;
    1108:           points[p++] = point;
    1109:         }
    1110:         const int cDim = this->section.getConstraintDimension(point);

    1112:         if (!cDim) {
    1113:           this->getUnconstrainedIndices(point, orientation);
    1114:         } else {
    1115:           this->getConstrainedIndices(point, orientation);
    1116:         }
    1117:       }
    1118:       template<typename Arrow>
    1119:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    1120:     public:
    1121:       const value_type *getValues() const {return this->values;};
    1122:       int  getSize() const {return this->i;};
    1123:       int  getMaxSize() const {return this->size;};
    1124:       void ensureSize(const int size) {
    1125:         this->clear();
    1126:         if (size > this->size) {
    1127:           this->size = size;
    1128:           if (this->allocated) {delete [] this->values;}
    1129:           this->values = new value_type[this->size];
    1130:           this->allocated = true;
    1131:           if (this->allocatedPoints) {delete [] this->points;}
    1132:           this->points = new point_type[this->size];
    1133:           this->allocatedPoints = true;
    1134:         }
    1135:       };
    1136:       void clear() {this->i = 0; this->p = 0;};
    1137:     };
    1138:     template<typename Order, typename Value>
    1139:     class IndicesVisitor<PetscSection, Order, Value> {
    1140:     public:
    1141:       typedef Value                      value_type;
    1142:       typedef typename Order::point_type point_type;
    1143:     protected:
    1144:       const PetscSection& section;
    1145:       // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
    1146:       Order&              order;
    1147:       int                 size;
    1148:       int                 i, p;
    1149:       bool                setBC;           // If true, returns indices for constrained dofs, otherwise negative values are returned
    1150:       //bool                skipConstraints; // If true, do not return constrained indices at all
    1151:       value_type         *values;
    1152:       bool                allocated;
    1153:       point_type         *points;
    1154:       PetscInt            nF;
    1155:       PetscInt           *fieldSize;
    1156:       PetscInt           *j;
    1157:     protected:
    1158:       void updatePoint(const point_type& point, const bool setBC, const int orientation = 1) {
    1159:         PetscInt        dim;  // The number of dof on this point
    1160:         PetscInt        cDim; // The nubmer of constraints on this point
    1161:         const PetscInt *cDof; // The indices of the constrained dofs on this point
    1162:         PetscInt        offset = this->order.getIndex(point);
    1163:         PetscInt        cInd   = 0;
    1164:         PetscErrorCode  ierr;

    1166:         PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
    1167:         PetscSectionGetConstraintDof(section, point, &cDim);CHKERRXX(ierr);
    1168:         if (!cDim || setBC) {
    1169:           if (orientation >= 0) {
    1170:             for(PetscInt k = 0; k < dim; ++k) {
    1171:               values[i+k] = offset+k;
    1172:             }
    1173:           } else {
    1174:             for(PetscInt k = 0; k < dim; ++k) {
    1175:               values[i+dim-k-1] = offset+k;
    1176:             }
    1177:           }
    1178:         } else {
    1179:           PetscSectionGetConstraintIndices(section, point, &cDof);CHKERRXX(ierr);
    1180:           if (orientation >= 0) {
    1181:             for(PetscInt k = 0; k < dim; ++k) {
    1182:               if ((cInd < cDim) && (k == cDof[cInd])) {
    1183:                 // Insert check for returning constrained indices
    1184:                 values[i+k] = -(offset+k+1);
    1185:                 ++cInd;
    1186:               } else {
    1187:                 values[i+k] = offset+k;
    1188:               }
    1189:             }
    1190:           } else {
    1191:             for(PetscInt k = 0; k < dim; ++k) {
    1192:               if ((cInd < cDim) && (k == cDof[cInd])) {
    1193:                 // Insert check for returning constrained indices
    1194:                 values[i+dim-k-1] = -(offset+k+1);
    1195:                 ++cInd;
    1196:               } else {
    1197:                 values[i+dim-k-1] = offset+k;
    1198:               }
    1199:             }
    1200:           }
    1201:         }
    1202:         i += dim;
    1203:       }
    1204:       void updatePointFields(const point_type& point, const bool setBC, const int orientation = 1) {
    1205:         PetscInt       offset = this->order.getIndex(point);
    1206:         PetscInt       fOff   = 0;

    1209:         for(PetscInt f = 0; f < nF; ++f) {
    1210:           PetscInt  dim;  // The number of dof for field f on this point
    1211:           PetscInt  comp; // The number of components for field f on this point
    1212:           PetscInt  cDim; // The nubmer of constraints for field f on this point
    1213:           const PetscInt *cDof; // The indices of the constrained dofs for field f on this point
    1214:           PetscInt  cInd = 0;

    1216:           PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
    1217:           PetscSectionGetFieldDof(section, point, f, &dim);CHKERRXX(ierr);
    1218:           PetscSectionGetFieldConstraintDof(section, point, f, &cDim);CHKERRXX(ierr);
    1219:           if (!cDim || setBC) {
    1220:             if (orientation >= 0) {
    1221:               for(PetscInt k = 0; k < dim; ++k) {
    1222:                 values[j[f]+k] = offset+fOff+k;
    1223:               }
    1224:             } else {
    1225:               for(PetscInt k = dim/comp-1; k >= 0; --k) {
    1226:                 for(PetscInt c = 0; c < comp; ++c) {
    1227:                   values[j[f]+(dim/comp-1-k)*comp+c] = offset+fOff+k*comp+c;
    1228:                 }
    1229:               }
    1230:             }
    1231:           } else {
    1232:             PetscSectionGetFieldConstraintIndices(section, point, f, &cDof);CHKERRXX(ierr);
    1233:             if (orientation >= 0) {
    1234:               for(PetscInt k = 0; k < dim; ++k) {
    1235:                 if ((cInd < cDim) && (k == cDof[cInd])) {
    1236:                   values[j[f]+k] = -(offset+fOff+k+1);
    1237:                   ++cInd;
    1238:                 } else {
    1239:                   values[j[f]+k] = offset+fOff+k;
    1240:                 }
    1241:               }
    1242:             } else {
    1243:               for(PetscInt k = dim/comp-1; k >= 0; --k) {
    1244:                 for(PetscInt c = 0; c < comp; ++c) {
    1245:                   PetscInt ind = k*comp+c;
    1246:                   if ((cInd < cDim) && (ind == cDof[cInd])) {
    1247:                     values[j[f]+(dim/comp-1-k)*comp+c] = -(offset+fOff+ind+1);
    1248:                     ++cInd;
    1249:                   } else {
    1250:                     values[j[f]+(dim/comp-1-k)*comp+c] = offset+fOff+ind;
    1251:                   }
    1252:                 }
    1253:               }
    1254:             }
    1255:           }
    1256:           fOff += dim - cDim;
    1257:           j[f] += dim;
    1258:           i    += dim;
    1259:         }
    1260:       }
    1261:     public:
    1262:       IndicesVisitor(const PetscSection& s, Order& o, const int size, const bool unique = false, const PetscInt fieldSize[] = NULL) : section(s), order(o), size(size), i(0), p(0), setBC(false) {

    1265:         PetscMalloc(this->size * sizeof(value_type), &this->values);CHKERRXX(ierr);
    1266:         this->allocated = true;
    1267:         this->points    = NULL;
    1268:         if (unique) {
    1269:           PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
    1270:         }
    1271:         nF = 0;
    1272:         this->fieldSize = this->j = NULL;
    1273:         if (fieldSize) {
    1274:           PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
    1275:           PetscMalloc2(nF,PetscInt,&this->fieldSize,nF,PetscInt,&j);CHKERRXX(ierr);
    1276:           for(PetscInt f = 0; f < nF; ++f) {
    1277:             this->fieldSize[f] = fieldSize[f];
    1278:           }
    1279:         }
    1280:         this->clear();
    1281:       };
    1282:       IndicesVisitor(const PetscSection& s, Order& o, const int size, value_type *values, const bool unique = false, const PetscInt fieldSize[] = NULL) : section(s), order(o), size(size), i(0), p(0), setBC(false) {

    1285:         this->values    = values;
    1286:         this->allocated = false;
    1287:         this->points    = NULL;
    1288:         if (unique) {
    1289:           PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
    1290:         }
    1291:         nF = 0;
    1292:         this->fieldSize = this->j = NULL;
    1293:         if (fieldSize) {
    1294:           PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
    1295:           PetscMalloc2(nF,PetscInt,&fieldSize,nF,PetscInt,&j);CHKERRXX(ierr);
    1296:           for(PetscInt f = 0; f < nF; ++f) {
    1297:             this->fieldSize[f] = fieldSize[f];
    1298:           }
    1299:         }
    1300:         this->clear();
    1301:       };
    1302:       ~IndicesVisitor() {
    1304:         if (this->allocated) {PetscFree(values);CHKERRXX(ierr);}
    1305:         PetscFree(points);CHKERRXX(ierr);
    1306:         PetscFree2(fieldSize,j);CHKERRXX(ierr);
    1307:       };
    1308:     public:
    1309:       inline void visitPoint(const point_type& point, const int orientation) {
    1310:         if (p >= size) {
    1311:           ostringstream msg;
    1312:           msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
    1313:           throw ALE::Exception(msg.str().c_str());
    1314:         }
    1315:         if (points) {
    1316:           PetscInt pp;
    1317:           for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
    1318:           if (pp != p) return;
    1319:           points[p++] = point;
    1320:         }
    1321:         if (nF) {
    1322:           updatePointFields(point, setBC, orientation);
    1323:         } else {
    1324:           updatePoint(point, setBC, orientation);
    1325:         }
    1326:       }
    1327:       template<typename Arrow>
    1328:       inline void visitArrow(const Arrow& arrow, const int orientation) {}
    1329:     public:
    1330:       const value_type *getValues() const {return this->values;};
    1331:       int  getSize() const {return this->i;};
    1332:       int  getMaxSize() const {return this->size;};
    1333:       void ensureSize(const int size) {
    1334:         this->clear();
    1335:         if (size > this->size) {

    1338:           this->size = size;
    1339:           if (this->allocated) {PetscFree(this->values);CHKERRXX(ierr);}
    1340:           PetscMalloc(this->size * sizeof(value_type), &this->values);CHKERRXX(ierr);
    1341:           this->allocated = true;
    1342:           PetscFree(this->points);CHKERRXX(ierr);
    1343:           PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
    1344:         }
    1345:       };
    1346:       void clear() {
    1347:         this->p = 0;
    1348:         this->i = 0;
    1349:         if (nF) {
    1350:           j[0] = 0;
    1351:           for(PetscInt f = 1; f < nF; ++f) {
    1352:             j[f] = j[f-1] + fieldSize[f-1];
    1353:           }
    1354:         }
    1355:       };
    1356:     };
    1357:     template<typename Sieve, typename Label>
    1358:     class MarkVisitor {
    1359:     protected:
    1360:       Label& label;
    1361:       int    marker;
    1362:     public:
    1363:       MarkVisitor(Label& l, const int marker) : label(l), marker(marker) {};
    1364:       inline void visitPoint(const typename Sieve::point_type& point) {
    1365:         this->label.setCone(this->marker, point);
    1366:       };
    1367:       inline void visitArrow(const typename Sieve::arrow_type&) {};
    1368:     };
    1369:   };

    1371:   template<typename Sieve>
    1372:   class ISieveTraversal {
    1373:   public:
    1374:     typedef typename Sieve::point_type point_type;
    1375:   public:
    1376:     template<typename Visitor>
    1377:     static void orientedClosure(const Sieve& sieve, const point_type& p, Visitor& v) {
    1378:       typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
    1379:       typedef ISieveVisitor::PointRetriever<Sieve,Retriever> TmpRetriever;
    1380:       Retriever    pV(200, v, true); // Correct estimate is pow(std::max(1, sieve->getMaxConeSize()), mesh->depth())
    1381:       TmpRetriever cV[2] = {TmpRetriever(200,pV), TmpRetriever(200,pV)};
    1382:       int          c     = 0;

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

    1388:       while(cV[c].getOrientedSize()) {
    1389:         const typename Retriever::oriented_point_type *cone     = cV[c].getOrientedPoints();
    1390:         const int                                      coneSize = cV[c].getOrientedSize();
    1391:         c = 1 - c;

    1393:         for(int p = 0; p < coneSize; ++p) {
    1394:           const typename Retriever::point_type& point     = cone[p].first;
    1395:           int                                   pO        = cone[p].second == 0 ? 1 : cone[p].second;
    1396:           const int                             pConeSize = sieve.getConeSize(point);

    1398:           if (pO < 0) {
    1399:             if (pO == -pConeSize) {
    1400:               sieve.orientedReverseCone(point, cV[c]);
    1401:             } else {
    1402:               const int numSkip = sieve.getConeSize(point) + pO;

    1404:               cV[c].setSkip(cV[c].getSize()+numSkip);
    1405:               cV[c].setLimit(cV[c].getSize()+pConeSize);
    1406:               sieve.orientedReverseCone(point, cV[c]);
    1407:               sieve.orientedReverseCone(point, cV[c]);
    1408:               cV[c].setSkip(0);
    1409:               cV[c].setLimit(0);
    1410:             }
    1411:           } else {
    1412:             if (pO == 1) {
    1413:               sieve.orientedCone(point, cV[c]);
    1414:             } else {
    1415:               const int numSkip = pO-1;

    1417:               cV[c].setSkip(cV[c].getSize()+numSkip);
    1418:               cV[c].setLimit(cV[c].getSize()+pConeSize);
    1419:               sieve.orientedCone(point, cV[c]);
    1420:               sieve.orientedCone(point, cV[c]);
    1421:               cV[c].setSkip(0);
    1422:               cV[c].setLimit(0);
    1423:             }
    1424:           }
    1425:         }
    1426:         cV[1-c].clear();
    1427:       }
    1428: #if 0
    1429:       // These contain arrows paired with orientations from the \emph{previous} arrow
    1430:       Obj<orientedArrowArray> cone    = new orientedArrowArray();
    1431:       Obj<orientedArrowArray> base    = new orientedArrowArray();
    1432:       coneSet                 seen;

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

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

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

    1449:           if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
    1450:             o = -(o+1);
    1451:           }
    1452:           if (o < 0) {
    1453:             const int size = pCone->size();

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

    1461:                   seen.insert(*c_iter);
    1462:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1463:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
    1464:                 }
    1465:               }
    1466:             } else {
    1467:               const int numSkip = size + o;
    1468:               int       count   = 0;

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

    1476:                   seen.insert(*c_iter);
    1477:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1478:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
    1479:                 }
    1480:               }
    1481:               count = 0;
    1482:               for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
    1483:                 if (count >= numSkip) break;
    1484:                 if (seen.find(*c_iter) == seen.end()) {
    1485:                   const arrow_type newArrow(*c_iter, arrow.source);
    1486:                   int              pointO = orientation->restrictPoint(newArrow)[0];

    1488:                   seen.insert(*c_iter);
    1489:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1490:                   closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
    1491:                 }
    1492:               }
    1493:             }
    1494:           } else {
    1495:             if (o == 1) {
    1496:               for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
    1497:                 if (seen.find(*c_iter) == seen.end()) {
    1498:                   const arrow_type newArrow(*c_iter, arrow.source);

    1500:                   seen.insert(*c_iter);
    1501:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1502:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
    1503:                 }
    1504:               }
    1505:             } else {
    1506:               const int numSkip = o-1;
    1507:               int       count   = 0;

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

    1514:                   seen.insert(*c_iter);
    1515:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1516:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
    1517:                 }
    1518:               }
    1519:               count = 0;
    1520:               for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
    1521:                 if (count >= numSkip) break;
    1522:                 if (seen.find(*c_iter) == seen.end()) {
    1523:                   const arrow_type newArrow(*c_iter, arrow.source);

    1525:                   seen.insert(*c_iter);
    1526:                   cone->push_back(oriented_arrow_type(newArrow, o));
    1527:                   closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
    1528:                 }
    1529:               }
    1530:             }
    1531:           }
    1532:         }
    1533:       }
    1534: #endif
    1535:     }
    1536:     template<typename Visitor>
    1537:     static void orientedStar(const Sieve& sieve, const point_type& p, Visitor& v) {
    1538:       typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
    1539:       Retriever sV[2] = {Retriever(200,v), Retriever(200,v)};
    1540:       int       s     = 0;

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

    1546:       while(sV[s].getOrientedSize()) {
    1547:         const typename Retriever::oriented_point_type *support     = sV[s].getOrientedPoints();
    1548:         const int                                      supportSize = sV[s].getOrientedSize();
    1549:         s = 1 - s;

    1551:         for(int p = 0; p < supportSize; ++p) {
    1552:           const typename Retriever::point_type& point        = support[p].first;
    1553:           int                                   pO           = support[p].second == 0 ? 1 : support[p].second;
    1554:           const int                             pSupportSize = sieve.getSupportSize(point);

    1556:           if (pO < 0) {
    1557:             if (pO == -pSupportSize) {
    1558:               sieve.orientedReverseSupport(point, sV[s]);
    1559:             } else {
    1560:               const int numSkip = sieve.getSupportSize(point) + pO;

    1562:               sV[s].setSkip(sV[s].getSize()+numSkip);
    1563:               sV[s].setLimit(sV[s].getSize()+pSupportSize);
    1564:               sieve.orientedReverseSupport(point, sV[s]);
    1565:               sieve.orientedReverseSupport(point, sV[s]);
    1566:               sV[s].setSkip(0);
    1567:               sV[s].setLimit(0);
    1568:             }
    1569:           } else {
    1570:             if (pO == 1) {
    1571:               sieve.orientedSupport(point, sV[s]);
    1572:             } else {
    1573:               const int numSkip = pO-1;

    1575:               sV[s].setSkip(sV[s].getSize()+numSkip);
    1576:               sV[s].setLimit(sV[s].getSize()+pSupportSize);
    1577:               sieve.orientedSupport(point, sV[s]);
    1578:               sieve.orientedSupport(point, sV[s]);
    1579:               sV[s].setSkip(0);
    1580:               sV[s].setLimit(0);
    1581:             }
    1582:           }
    1583:         }
    1584:         sV[1-s].clear();
    1585:       }
    1586:     }
    1587:   };

    1589:   namespace IFSieveDef {
    1590:     template<typename PointType_>
    1591:     class Sequence {
    1592:     public:
    1593:       typedef PointType_ point_type;
    1594:       class const_iterator {
    1595:       public:
    1596:         // Standard iterator typedefs
    1597:         typedef std::input_iterator_tag iterator_category;
    1598:         typedef PointType_              value_type;
    1599:         typedef int                     difference_type;
    1600:         typedef value_type*             pointer;
    1601:         typedef value_type&             reference;
    1602:       protected:
    1603:         const point_type *_data;
    1604:         int               _pos;
    1605:       public:
    1606:         const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
    1607:         virtual ~const_iterator() {};
    1608:       public:
    1609:         virtual bool              operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
    1610:         virtual bool              operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
    1611:         virtual const value_type  operator*() const {return this->_data[this->_pos];};
    1612:         virtual const_iterator&   operator++() {++this->_pos; return *this;};
    1613:         virtual const_iterator    operator++(int n) {
    1614:           const_iterator tmp(*this);
    1615:           ++this->_pos;
    1616:           return tmp;
    1617:         };
    1618:       };
    1619:       typedef const_iterator iterator;
    1620:     protected:
    1621:       const point_type *_data;
    1622:       int               _begin;
    1623:       int               _end;
    1624:     public:
    1625:       Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
    1626:       virtual ~Sequence() {};
    1627:     public:
    1628:       virtual iterator begin() const {return iterator(_data, _begin);};
    1629:       virtual iterator end()   const {return iterator(_data, _end);};
    1630:       size_t size() const {return(_end - _begin);}
    1631:       void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
    1632:     };
    1633:   }

    1635:   // Interval Final Sieve
    1636:   // This is just two CSR matrices that give cones and supports
    1637:   //   It is completely static and cannot be resized
    1638:   //   It will operator on visitors, rather than sequences (which are messy)
    1639:   template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
    1640:   class IFSieve : public ParallelObject {
    1641:   public:
    1642:     // Types
    1643:     typedef IFSieve<Point_,Allocator_>         this_type;
    1644:     typedef Point_                             point_type;
    1645:     typedef SimpleArrow<point_type,point_type> arrow_type;
    1646:     typedef typename arrow_type::source_type   source_type;
    1647:     typedef typename arrow_type::target_type   target_type;
    1648:     typedef int                                index_type;
    1649:     // Allocators
    1650:     typedef Allocator_                                                        point_allocator_type;
    1651:     typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
    1652:     typedef typename point_allocator_type::template rebind<int>::other        int_allocator_type;
    1653:     // Interval
    1654:     typedef Interval<point_type, point_allocator_type> chart_type;
    1655:     // Dynamic structure
    1656:     typedef std::map<point_type, std::vector<point_type> > newpoints_type;
    1657:     // Iterator interface
    1658:     typedef typename IFSieveDef::Sequence<point_type> coneSequence;
    1659:     typedef typename IFSieveDef::Sequence<point_type> supportSequence;
    1660:     // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
    1661:     typedef std::set<point_type>   pointSet;
    1662:     typedef ALE::array<point_type> pointArray;
    1663:     typedef pointSet               coneSet;
    1664:     typedef pointSet               supportSet;
    1665:     typedef pointArray             coneArray;
    1666:     typedef pointArray             supportArray;
    1667:   protected:
    1668:     // Arrow Containers
    1669:     typedef index_type *offsets_type;
    1670:     typedef point_type *cones_type;
    1671:     typedef point_type *supports_type;
    1672:     // Decorators
    1673:     typedef int        *orientations_type;
    1674:   protected:
    1675:     // Data
    1676:     bool                 indexAllocated;
    1677:     offsets_type         coneOffsets;
    1678:     offsets_type         supportOffsets;
    1679:     bool                 pointAllocated;
    1680:     index_type           maxConeSize;
    1681:     index_type           maxSupportSize;
    1682:     index_type           baseSize;
    1683:     index_type           capSize;
    1684:     cones_type           cones;
    1685:     supports_type        supports;
    1686:     bool                 orientCones;
    1687:     orientations_type    coneOrientations;
    1688:     chart_type           chart;
    1689:     int_allocator_type   intAlloc;
    1690:     index_allocator_type indexAlloc;
    1691:     point_allocator_type pointAlloc;
    1692:     newpoints_type       newCones;
    1693:     newpoints_type       newSupports;
    1694:     // Sequences
    1695:     Obj<coneSequence>    coneSeq;
    1696:     Obj<supportSequence> supportSeq;
    1697:   protected: // Memory Management
    1698:     void createIndices() {
    1699:       this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
    1700:       this->coneOffsets -= this->chart.min();
    1701:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
    1702:       this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
    1703:       this->supportOffsets -= this->chart.min();
    1704:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
    1705:       this->indexAllocated = true;
    1706:     };
    1707:     void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
    1708:       if (*coneOffsets) {
    1709:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
    1710:         *coneOffsets += chart.min();
    1711:         indexAlloc.deallocate(*coneOffsets, chart.size()+1);
    1712:         *coneOffsets = NULL;
    1713:       }
    1714:       if (*supportOffsets) {
    1715:         for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
    1716:         *supportOffsets += chart.min();
    1717:         indexAlloc.deallocate(*supportOffsets, chart.size()+1);
    1718:         *supportOffsets = NULL;
    1719:       }
    1720:     };
    1721:     void destroyIndices() {
    1722:       this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
    1723:       this->indexAllocated = false;
    1724:       this->maxConeSize    = 0;
    1725:       this->maxSupportSize = 0;
    1726:       this->baseSize       = -1;
    1727:       this->capSize        = -1;
    1728:     };
    1729:     void createPoints() {
    1730:       this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
    1731:       for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
    1732:       this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
    1733:       for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
    1734:       if (orientCones) {
    1735:         this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
    1736:         for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
    1737:       }
    1738:       this->pointAllocated = true;
    1739:     };
    1740:     void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
    1741:       if (*cones) {
    1742:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
    1743:         pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
    1744:         *cones = NULL;
    1745:       }
    1746:       if (*supports) {
    1747:         for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
    1748:         pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
    1749:         *supports = NULL;
    1750:       }
    1751:       if (*coneOrientations) {
    1752:         for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
    1753:         intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
    1754:         *coneOrientations = NULL;
    1755:       }
    1756:     };
    1757:     void destroyPoints() {
    1758:       this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
    1759:       this->pointAllocated = false;
    1760:     };
    1761:     void prefixSum(const offsets_type array) {
    1762:       for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
    1763:         array[p] = array[p] + array[p-1];
    1764:       }
    1765:     };
    1766:     void calculateBaseAndCapSize() {
    1767:       this->baseSize = 0;
    1768:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    1769:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
    1770:           ++this->baseSize;
    1771:         }
    1772:       }
    1773:       this->capSize = 0;
    1774:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    1775:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
    1776:           ++this->capSize;
    1777:         }
    1778:       }
    1779:     };
    1780:   public:
    1781:     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) {
    1782:       this->coneSeq    = new coneSequence(NULL, 0, 0);
    1783:       this->supportSeq = new supportSequence(NULL, 0, 0);
    1784:     };
    1785:     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) {
    1786:       this->setChart(chart_type(min, max));
    1787:       this->coneSeq    = new coneSequence(NULL, 0, 0);
    1788:       this->supportSeq = new supportSequence(NULL, 0, 0);
    1789:     };
    1790:     ~IFSieve() {
    1791:       this->destroyPoints();
    1792:       this->destroyIndices();
    1793:     };
    1794:   public: // Accessors
    1795:     const chart_type& getChart() const {return this->chart;};
    1796:     void setChart(const chart_type& chart) {
    1797:       this->destroyPoints();
    1798:       this->destroyIndices();
    1799:       this->chart = chart;
    1800:       this->createIndices();
    1801:     };
    1802:     index_type getMaxConeSize() const {
    1803:       return this->maxConeSize;
    1804:     };
    1805:     index_type getMaxSupportSize() const {
    1806:       return this->maxSupportSize;
    1807:     };
    1808:     bool baseContains(const point_type& p) const {
    1809:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    1810:       this->chart.checkPoint(p);

    1812:       if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
    1813:         return true;
    1814:       }
    1815:       return false;
    1816:     };
    1817:     bool orientedCones() const {return this->orientCones;};
    1818:     // Raw array access
    1819:     offsets_type      getConeOffsets() {return this->coneOffsets;};
    1820:     offsets_type      getSupportOffsets() {return this->supportOffsets;};
    1821:     cones_type        getCones() {return this->cones;};
    1822:     supports_type     getSupports() {return this->supports;};
    1823:     orientations_type getConeOrientations() {return this->coneOrientations;};
    1824:   public: // Construction
    1825:     index_type getConeSize(const point_type& p) const {
    1826:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    1827:       this->chart.checkPoint(p);
    1828:       return this->coneOffsets[p+1]-this->coneOffsets[p];
    1829:     };
    1830:     void setConeSize(const point_type& p, const index_type c) {
    1831:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
    1832:       this->chart.checkPoint(p);
    1833:       this->coneOffsets[p+1] = c;
    1834:       this->maxConeSize = std::max(this->maxConeSize, c);
    1835:     };
    1836:     void addConeSize(const point_type& p, const index_type c) {
    1837:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
    1838:       this->chart.checkPoint(p);
    1839:       this->coneOffsets[p+1] += c;
    1840:       this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
    1841:     };
    1842:     index_type getSupportSize(const point_type& p) const {
    1843:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    1844:       this->chart.checkPoint(p);
    1845:       return this->supportOffsets[p+1]-this->supportOffsets[p];
    1846:     };
    1847:     void setSupportSize(const point_type& p, const index_type s) {
    1848:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
    1849:       this->chart.checkPoint(p);
    1850:       this->supportOffsets[p+1] = s;
    1851:       this->maxSupportSize = std::max(this->maxSupportSize, s);
    1852:     };
    1853:     void addSupportSize(const point_type& p, const index_type s) {
    1854:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
    1855:       this->chart.checkPoint(p);
    1856:       this->supportOffsets[p+1] += s;
    1857:       this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
    1858:     };
    1859:     void allocate() {
    1860:       if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
    1861:       this->prefixSum(this->coneOffsets);
    1862:       this->prefixSum(this->supportOffsets);
    1863:       this->createPoints();
    1864:       this->calculateBaseAndCapSize();
    1865:     }
    1866:     // Purely for backwards compatibility
    1867:     template<typename Color>
    1868:     void addArrow(const point_type& p, const point_type& q, const Color c, const bool forceSupport = false) {
    1869:       this->addArrow(p, q, forceSupport);
    1870:     }
    1871:     void addArrow(const point_type& p, const point_type& q, const bool forceSupport = false) {
    1872:       if (!this->chart.hasPoint(q)) {
    1873:         if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
    1874:           const index_type start = this->coneOffsets[q];
    1875:           const index_type end   = this->coneOffsets[q+1];

    1877:           for(int c = start; c < end; ++c) {
    1878:             this->newCones[q].push_back(this->cones[c]);
    1879:           }
    1880:         }
    1881:         this->newCones[q].push_back(p);
    1882:       }
    1883:       if (!this->chart.hasPoint(p) || forceSupport) {
    1884:         if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
    1885:           const index_type start = this->supportOffsets[p];
    1886:           const index_type end   = this->supportOffsets[p+1];

    1888:           for(int s = start; s < end; ++s) {
    1889:             this->newSupports[p].push_back(this->supports[s]);
    1890:           }
    1891:         }
    1892:         this->newSupports[p].push_back(q);
    1893:       }
    1894:     };
    1895:     void reallocate() {
    1896:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    1897:       if (!this->newCones.size() && !this->newSupports.size()) return;
    1898:       const chart_type     oldChart            = this->chart;
    1899:       offsets_type         oldConeOffsets      = this->coneOffsets;
    1900:       offsets_type         oldSupportOffsets   = this->supportOffsets;
    1901:       cones_type           oldCones            = this->cones;
    1902:       supports_type        oldSupports         = this->supports;
    1903:       orientations_type    oldConeOrientations = this->coneOrientations;
    1904:       point_type           min                 = this->chart.min();
    1905:       point_type           max                 = this->chart.max()-1;

    1907:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
    1908:         min = std::min(min, c_iter->first);
    1909:         max = std::max(max, c_iter->first);
    1910:       }
    1911:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
    1912:         min = std::min(min, s_iter->first);
    1913:         max = std::max(max, s_iter->first);
    1914:       }
    1915:       this->chart = chart_type(min, max+1);
    1916:       this->createIndices();
    1917:       // Copy sizes (converted from offsets)
    1918:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
    1919:         this->coneOffsets[p+1]    = oldConeOffsets[p+1]-oldConeOffsets[p];
    1920:         this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
    1921:       }
    1922:       // Inject new sizes
    1923:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
    1924:         this->coneOffsets[c_iter->first+1]    = c_iter->second.size();
    1925:         this->maxConeSize                     = std::max(this->maxConeSize,    (int) c_iter->second.size());
    1926:       }
    1927:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
    1928:         this->supportOffsets[s_iter->first+1] = s_iter->second.size();
    1929:         this->maxSupportSize                  = std::max(this->maxSupportSize, (int) s_iter->second.size());
    1930:       }
    1931:       this->prefixSum(this->coneOffsets);
    1932:       this->prefixSum(this->supportOffsets);
    1933:       this->createPoints();
    1934:       this->calculateBaseAndCapSize();
    1935:       // Copy cones and supports
    1936:       for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
    1937:         const index_type cStart  = this->coneOffsets[p];
    1938:         const index_type cOStart = oldConeOffsets[p];
    1939:         const index_type cOEnd   = oldConeOffsets[p+1];
    1940:         const index_type sStart  = this->supportOffsets[p];
    1941:         const index_type sOStart = oldSupportOffsets[p];
    1942:         const index_type sOEnd   = oldSupportOffsets[p+1];

    1944:         for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
    1945:           this->cones[c] = oldCones[cO];
    1946:         }
    1947:         for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
    1948:           this->supports[s] = oldSupports[sO];
    1949:         }
    1950:         if (this->orientCones) {
    1951:           for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
    1952:             this->coneOrientations[c] = oldConeOrientations[cO];
    1953:           }
    1954:         }
    1955:       }
    1956:       // Inject new cones and supports
    1957:       for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
    1958:         index_type start = this->coneOffsets[c_iter->first];

    1960:         for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
    1961:           this->cones[start++] = *p_iter;
    1962:         }
    1963:         if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
    1964:       }
    1965:       for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
    1966:         index_type start = this->supportOffsets[s_iter->first];

    1968:         for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
    1969:           this->supports[start++] = *p_iter;
    1970:         }
    1971:         if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
    1972:       }
    1973:       this->newCones.clear();
    1974:       this->newSupports.clear();
    1975:       this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
    1976:       this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
    1977:     };
    1978:     // Recalculate the support offsets and fill the supports
    1979:     //   This is used when an IFSieve is being used as a label
    1980:     void recalculateLabel() {
    1981:       ISieveVisitor::PointRetriever<IFSieve> v(1);

    1983:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
    1984:         this->supportOffsets[p+1] = 0;
    1985:       }
    1986:       this->maxSupportSize = 0;
    1987:       for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
    1988:         this->cone(p, v);
    1989:         if (v.getSize()) {
    1990:           this->supportOffsets[v.getPoints()[0]+1]++;
    1991:           this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
    1992:         }
    1993:         v.clear();
    1994:       }
    1995:       this->prefixSum(this->supportOffsets);
    1996:       this->calculateBaseAndCapSize();
    1997:       this->symmetrize();
    1998:     };
    1999:     void setCone(const point_type cone[], const point_type& p) {
    2000:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2001:       this->chart.checkPoint(p);
    2002:       const index_type start = this->coneOffsets[p];
    2003:       const index_type end   = this->coneOffsets[p+1];

    2005:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
    2006:         this->cones[c] = cone[i];
    2007:       }
    2008:     };
    2009:     void setCone(const point_type cone, const point_type& p) {
    2010:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2011:       this->chart.checkPoint(p);
    2012:       const index_type start = this->coneOffsets[p];
    2013:       const index_type end   = this->coneOffsets[p+1];

    2015:       if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
    2016:       this->cones[start] = cone;
    2017:     };
    2018: #if 0
    2019:     template<typename PointSequence>
    2020:     void setCone(const PointSequence& cone, const point_type& p) {
    2021:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2022:       this->chart.checkPoint(p);
    2023:       const index_type start = this->coneOffsets[p];
    2024:       const index_type end   = this->coneOffsets[p+1];
    2025:       if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
    2026:       typename PointSequence::iterator c_iter = cone.begin();

    2028:       for(index_type c = start; c < end; ++c, ++c_iter) {
    2029:         this->cones[c] = *c_iter;
    2030:       }
    2031:     };
    2032: #endif
    2033:     void setSupport(const point_type& p, const point_type support[]) {
    2034:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2035:       this->chart.checkPoint(p);
    2036:       const index_type start = this->supportOffsets[p];
    2037:       const index_type end   = this->supportOffsets[p+1];

    2039:       for(index_type s = start, i = 0; s < end; ++s, ++i) {
    2040:         this->supports[s] = support[i];
    2041:       }
    2042:     };
    2043: #if 0
    2044:     template<typename PointSequence>
    2045:     void setSupport(const point_type& p, const PointSequence& support) {
    2046:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2047:       this->chart.checkPoint(p);
    2048:       const index_type start = this->supportOffsets[p];
    2049:       const index_type end   = this->supportOffsets[p+1];
    2050:       if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
    2051:       typename PointSequence::iterator s_iter = support.begin();

    2053:       for(index_type s = start; s < end; ++s, ++s_iter) {
    2054:         this->supports[s] = *s_iter;
    2055:       }
    2056:     };
    2057: #endif
    2058:     void setConeOrientation(const int coneOrientation[], const point_type& p) {
    2059:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2060:       this->chart.checkPoint(p);
    2061:       const index_type start = this->coneOffsets[p];
    2062:       const index_type end   = this->coneOffsets[p+1];

    2064:       for(index_type c = start, i = 0; c < end; ++c, ++i) {
    2065:         this->coneOrientations[c] = coneOrientation[i];
    2066:       }
    2067:     };
    2068:     void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
    2069:       for(point_type p = 0; p < numCells; ++p) {
    2070:         const index_type start = p*numCorners;
    2071:         const index_type end   = (p+1)*numCorners;

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

    2076:           this->supportOffsets[q+1]++;
    2077:         }
    2078:       }
    2079:       for(point_type p = numCells; p < this->chart.max(); ++p) {
    2080:         this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
    2081:       }
    2082:     };
    2083:     void symmetrize() {
    2084:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
    2085:       offsets -= this->chart.min();
    2086:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
    2087:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

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

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

    2096:           this->chart.checkPoint(q);
    2097:           this->supports[this->supportOffsets[q]+offsets[q]] = p;
    2098:           ++offsets[q];
    2099:         }
    2100:       }
    2101:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
    2102:       indexAlloc.deallocate(offsets, this->chart.size()+1);
    2103:     }
    2104:     index_type getBaseSize() const {
    2105:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2106:       return this->baseSize;
    2107:     }
    2108:     index_type getCapSize() const {
    2109:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2110:       return this->capSize;
    2111:     }
    2112:     template<typename _Section>
    2113:     void relabel(_Section& labeling) {

    2116:       index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
    2117:       offsets -= this->chart.min();
    2118:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
    2119:       // Recalculate coneOffsets
    2120:       for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2121:         const point_type newP = labeling.restrictPoint(p)[0];

    2123:         offsets[newP+1] = this->getConeSize(p);
    2124:       }
    2125:       this->prefixSum(offsets);
    2126:       PetscMemcpy(this->coneOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
    2127:       // Recalculate supportOffsets
    2128:       for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2129:         const point_type newP = labeling.restrictPoint(p)[0];

    2131:         offsets[newP+1] = this->getSupportSize(p);
    2132:       }
    2133:       this->prefixSum(offsets);
    2134:       PetscMemcpy(this->supportOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
    2135:       for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
    2136:       indexAlloc.deallocate(offsets, this->chart.size()+1);
    2137:       index_type  size = std::max(this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()],
    2138:                                   this->supportOffsets[this->chart.max()] - this->supportOffsets[this->chart.min()]);
    2139:       index_type *orientations = offsets = indexAlloc.allocate(size);
    2140:       for(index_type i = 0; i < size; ++i) {indexAlloc.construct(orientations+i, index_type(0));}
    2141:       // Recalculate coneOrientations
    2142:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
    2143:         const point_type newP  = labeling.restrictPoint(p)[0];
    2144:         const index_type start = this->coneOffsets[newP];
    2145:         const index_type end   = this->coneOffsets[newP+1];

    2147:         for(index_type c = start; c < end; ++c, ++offset) {
    2148:           orientations[c] = this->coneOrientations[offset];
    2149:         }
    2150:       }
    2151:       PetscMemcpy(this->coneOrientations, orientations, (this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()])*sizeof(index_type));CHKERRXX(ierr);
    2152:       for(index_type i = 0; i < size; ++i) {indexAlloc.destroy(orientations+i);}
    2153:       indexAlloc.deallocate(orientations, size);
    2154:       // Recalculate cones
    2155:       point_type *array = pointAlloc.allocate(size);

    2157:       for(index_type i = 0; i < size; ++i) {pointAlloc.construct(array+i, point_type(0));}
    2158:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
    2159:         const point_type newP  = labeling.restrictPoint(p)[0];
    2160:         const index_type start = this->coneOffsets[newP];
    2161:         const index_type end   = this->coneOffsets[newP+1];

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

    2166:           array[c] = newQ;
    2167:         }
    2168:       }
    2169:       PetscMemcpy(this->cones, array, size*sizeof(point_type));CHKERRXX(ierr);
    2170:       // Recalculate supports
    2171:       for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
    2172:         const point_type newP  = labeling.restrictPoint(p)[0];
    2173:         const index_type start = this->supportOffsets[newP];
    2174:         const index_type end   = this->supportOffsets[newP+1];

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

    2179:           array[c] = newQ;
    2180:         }
    2181:       }
    2182:       PetscMemcpy(this->supports, array, size*sizeof(point_type));CHKERRXX(ierr);
    2183:       for(index_type i = 0; i < size; ++i) {pointAlloc.destroy(array+i);}
    2184:       pointAlloc.deallocate(array, size);
    2185:     }
    2186:   public: // Traversals
    2187:     template<typename Visitor>
    2188:     void roots(const Visitor& v) const {
    2189:       this->roots(const_cast<Visitor&>(v));
    2190:     }
    2191:     template<typename Visitor>
    2192:     void roots(Visitor& v) const {
    2193:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

    2195:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2196:         if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
    2197:           if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
    2198:             v.visitPoint(p);
    2199:           }
    2200:         }
    2201:       }
    2202:     }
    2203:     int numRoots() {
    2204:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2205:       int n = 0;

    2207:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2208:         if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
    2209:           if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
    2210:             ++n;
    2211:           }
    2212:         }
    2213:       }
    2214:       return n;
    2215:     }
    2216:     template<typename Visitor>
    2217:     void leaves(const Visitor& v) const {
    2218:       this->leaves(const_cast<Visitor&>(v));
    2219:     }
    2220:     template<typename Visitor>
    2221:     void leaves(Visitor& v) const {
    2222:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

    2224:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2225:         if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
    2226:           if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
    2227:             v.visitPoint(p);
    2228:           }
    2229:         }
    2230:       }
    2231:     }
    2232:     int numLeaves() {
    2233:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2234:       int n = 0;

    2236:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2237:         if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
    2238:           if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
    2239:             ++n;
    2240:           }
    2241:         }
    2242:       }
    2243:       return n;
    2244:     }
    2245:     template<typename Visitor>
    2246:     void base(const Visitor& v) const {
    2247:       this->base(const_cast<Visitor&>(v));
    2248:     }
    2249:     template<typename Visitor>
    2250:     void base(Visitor& v) const {
    2251:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

    2253:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2254:         if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
    2255:           v.visitPoint(p);
    2256:         }
    2257:       }
    2258:     }
    2259:     template<typename Visitor>
    2260:     void cap(const Visitor& v) const {
    2261:       this->cap(const_cast<Visitor&>(v));
    2262:     }
    2263:     template<typename Visitor>
    2264:     void cap(Visitor& v) const {
    2265:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}

    2267:       for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
    2268:         if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
    2269:           v.visitPoint(p);
    2270:         }
    2271:       }
    2272:     }
    2273:     template<typename PointSequence, typename Visitor>
    2274:     void cone(const PointSequence& points, Visitor& v) const {
    2275:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2276:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    2277:         const point_type p = *p_iter;
    2278:         this->chart.checkPoint(p);
    2279:         const index_type start = this->coneOffsets[p];
    2280:         const index_type end   = this->coneOffsets[p+1];

    2282:         for(index_type c = start; c < end; ++c) {
    2283:           v.visitArrow(arrow_type(this->cones[c], p));
    2284:           v.visitPoint(this->cones[c]);
    2285:         }
    2286:       }
    2287:     }
    2288:     template<typename Visitor>
    2289:     void cone(const point_type& p, Visitor& v) const {
    2290:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2291:       this->chart.checkPoint(p);
    2292:       const index_type start = this->coneOffsets[p];
    2293:       const index_type end   = this->coneOffsets[p+1];

    2295:       for(index_type c = start; c < end; ++c) {
    2296:         v.visitArrow(arrow_type(this->cones[c], p));
    2297:         v.visitPoint(this->cones[c]);
    2298:       }
    2299:     }
    2300:     const Obj<coneSequence>& cone(const point_type& p) const {
    2301:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2302:       if (!this->chart.hasPoint(p)) {
    2303:         this->coneSeq->reset(this->cones, 0, 0);
    2304:       } else {
    2305:         this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
    2306:       }
    2307:       return this->coneSeq;
    2308:     }
    2309:     template<typename PointSequence, typename Visitor>
    2310:     void support(const PointSequence& points, Visitor& v) const {
    2311:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2312:       for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    2313:         const point_type p = *p_iter;
    2314:         this->chart.checkPoint(p);
    2315:         const index_type start = this->supportOffsets[p];
    2316:         const index_type end   = this->supportOffsets[p+1];

    2318:         for(index_type s = start; s < end; ++s) {
    2319:           v.visitArrow(arrow_type(p, this->supports[s]));
    2320:           v.visitPoint(this->supports[s]);
    2321:         }
    2322:       }
    2323:     }
    2324:     template<typename Visitor>
    2325:     void support(const point_type& p, Visitor& v) const {
    2326:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2327:       this->chart.checkPoint(p);
    2328:       const index_type start = this->supportOffsets[p];
    2329:       const index_type end   = this->supportOffsets[p+1];

    2331:       for(index_type s = start; s < end; ++s) {
    2332:         v.visitArrow(arrow_type(p, this->supports[s]));
    2333:         v.visitPoint(this->supports[s]);
    2334:       }
    2335:     }
    2336:     const Obj<supportSequence>& support(const point_type& p) const {
    2337:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2338:       if (!this->chart.hasPoint(p)) {
    2339:         this->supportSeq->reset(this->supports, 0, 0);
    2340:       } else {
    2341:         this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
    2342:       }
    2343:       return this->supportSeq;
    2344:     }
    2345:     template<typename Visitor>
    2346:     void orientedCone(const point_type& p, Visitor& v) const {
    2347:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
    2348:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2349:       this->chart.checkPoint(p);
    2350:       const index_type start = this->coneOffsets[p];
    2351:       const index_type end   = this->coneOffsets[p+1];

    2353:       for(index_type c = start; c < end; ++c) {
    2354:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
    2355:         v.visitPoint(this->cones[c], this->coneOrientations[c]);
    2356:       }
    2357:     }
    2358:     template<typename Visitor>
    2359:     void orientedReverseCone(const point_type& p, Visitor& v) const {
    2360:       if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
    2361:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2362:       this->chart.checkPoint(p);
    2363:       const index_type start = this->coneOffsets[p];
    2364:       const index_type end   = this->coneOffsets[p+1];

    2366:       for(index_type c = end-1; c >= start; --c) {
    2367:         v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
    2368:         v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
    2369:       }
    2370:     }
    2371:     template<typename Visitor>
    2372:     void orientedSupport(const point_type& p, Visitor& v) const {
    2373:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
    2374:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2375:       this->chart.checkPoint(p);
    2376:       const index_type start = this->supportOffsets[p];
    2377:       const index_type end   = this->supportOffsets[p+1];

    2379:       for(index_type s = start; s < end; ++s) {
    2380:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
    2381:         //v.visitPoint(this->supports[s], this->supportOrientations[s]);
    2382:         v.visitArrow(arrow_type(this->supports[s], p), 0);
    2383:         v.visitPoint(this->supports[s], 0);
    2384:       }
    2385:     }
    2386:     template<typename Visitor>
    2387:     void orientedReverseSupport(const point_type& p, Visitor& v) const {
    2388:       //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
    2389:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2390:       this->chart.checkPoint(p);
    2391:       const index_type start = this->supportOffsets[p];
    2392:       const index_type end   = this->supportOffsets[p+1];

    2394:       for(index_type s = end-1; s >= start; --s) {
    2395:         //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
    2396:         //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
    2397:         v.visitArrow(arrow_type(this->supports[s], p), 0);
    2398:         v.visitPoint(this->supports[s], 0);
    2399:       }
    2400:     }
    2401:     // Currently does only 1 level
    2402:     //   Does not check for uniqueness
    2403:     template<typename Visitor>
    2404:     void meet(const point_type& p, const point_type& q, Visitor& v) const {
    2405:       if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
    2406:       this->chart.checkPoint(p);
    2407:       this->chart.checkPoint(q);
    2408:       const index_type startP = this->coneOffsets[p];
    2409:       const index_type endP   = this->coneOffsets[p+1];
    2410:       const index_type startQ = this->coneOffsets[q];
    2411:       const index_type endQ   = this->coneOffsets[q+1];

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

    2416:         for(index_type cQ = startQ; cQ < endQ; ++cQ) {
    2417:           if (c1 == this->cones[cQ]) v.visitPoint(c1);
    2418:         }
    2419:         if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
    2420:       }
    2421:     }
    2422:     // Currently does only 1 level
    2423:     template<typename Sequence, typename Visitor>
    2424:     void join(const Sequence& points, Visitor& v) const {
    2425:       typedef std::set<point_type> pointSet;
    2426:       pointSet intersect[2] = {pointSet(), pointSet()};
    2427:       pointSet tmp;
    2428:       int      p = 0;
    2429:       int      c = 0;

    2431:       for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    2432:         this->chart.checkPoint(*p_iter);
    2433:         tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
    2434:         if (p == 0) {
    2435:           intersect[1-c].insert(tmp.begin(), tmp.end());
    2436:           p++;
    2437:         } else {
    2438:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
    2439:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
    2440:           intersect[c].clear();
    2441:         }
    2442:         c = 1 - c;
    2443:         tmp.clear();
    2444:       }
    2445:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
    2446:         v.visitPoint(*p_iter);
    2447:       }
    2448:     }
    2449:     // Helper function
    2450:     void insertNSupport(point_type p, pointSet& set, const int depth) {
    2451:       const index_type start = this->supportOffsets[p];
    2452:       const index_type end   = this->supportOffsets[p+1];

    2454:       if (depth == 1) {
    2455:         set.insert(&this->supports[start], &this->supports[end]);
    2456:       } else {
    2457:         for(index_type s = start; s < end; ++s) {
    2458:           this->insertNSupport(this->supports[s], set, depth-1);
    2459:         }
    2460:       }
    2461:     }
    2462:     // Gives only the join of depth n
    2463:     template<typename SequenceIterator, typename Visitor>
    2464:     void nJoin(const SequenceIterator& pointsBegin, const SequenceIterator& pointsEnd, const int depth, Visitor& v) {
    2465:       typedef std::set<point_type> pointSet;
    2466:       pointSet intersect[2] = {pointSet(), pointSet()};
    2467:       pointSet tmp;
    2468:       int      p = 0;
    2469:       int      c = 0;

    2471:       for(SequenceIterator p_iter = pointsBegin; p_iter != pointsEnd; ++p_iter) {
    2472:         this->chart.checkPoint(*p_iter);
    2473:         // Put points in the nSupport into tmp (duplicates are fine since it is a set)
    2474:         this->insertNSupport(*p_iter, tmp, depth);
    2475:         if (p == 0) {
    2476:           intersect[1-c].insert(tmp.begin(), tmp.end());
    2477:           p++;
    2478:         } else {
    2479:           std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
    2480:                                 std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
    2481:           intersect[c].clear();
    2482:         }
    2483:         c = 1 - c;
    2484:         tmp.clear();
    2485:       }
    2486:       for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
    2487:         v.visitPoint(*p_iter);
    2488:       }
    2489:     }
    2490:   public: // Viewing
    2491:     void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
    2492:       ostringstream txt;
    2493:       int rank;

    2495:       if (comm == MPI_COMM_NULL) {
    2496:         comm = this->comm();
    2497:         rank = this->commRank();
    2498:       } else {
    2499:         MPI_Comm_rank(comm, &rank);
    2500:       }
    2501:       if (name == "") {
    2502:         if(rank == 0) {
    2503:           txt << "viewing an IFSieve" << std::endl;
    2504:         }
    2505:       } else {
    2506:         if(rank == 0) {
    2507:           txt << "viewing IFSieve '" << name << "'" << std::endl;
    2508:         }
    2509:       }
    2510:       PetscSynchronizedPrintf(comm, "Max sizes cone: %d support: %d\n", this->getMaxConeSize(), this->getMaxSupportSize());
    2511:       if(rank == 0) {
    2512:         txt << "cap --> base:" << std::endl;
    2513:       }
    2514:       ISieveVisitor::PrintVisitor pV(txt, rank);
    2515:       this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
    2516:       PetscSynchronizedPrintf(comm, txt.str().c_str());
    2517:       PetscSynchronizedFlush(comm);
    2518:       ostringstream txt2;

    2520:       if(rank == 0) {
    2521:         txt2 << "base <-- cap:" << std::endl;
    2522:       }
    2523:       ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
    2524:       this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
    2525:       PetscSynchronizedPrintf(comm, txt2.str().c_str());
    2526:       PetscSynchronizedFlush(comm);
    2527:       if (orientCones) {
    2528:         ostringstream txt3;

    2530:         if(rank == 0) {
    2531:           txt3 << "Orientation:" << std::endl;
    2532:         }
    2533:         ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
    2534:         this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
    2535:         PetscSynchronizedPrintf(comm, txt3.str().c_str());
    2536:         PetscSynchronizedFlush(comm);
    2537:       }
    2538:     };
    2539:   };

    2541:   class ISieveConverter {
    2542:   public:
    2543:     template<typename Sieve, typename ISieve, typename Renumbering>
    2544:     static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
    2545:       // First construct a renumbering of the sieve points
    2546:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
    2547:       const Obj<typename Sieve::capSequence>&  cap  = sieve.cap();
    2548:       typename ISieve::point_type              min  = 0;
    2549:       typename ISieve::point_type              max  = 0;

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

    2554:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2555:           if (sieve.support(*b_iter)->size() == 0) {
    2556:             renumbering[*b_iter] = max++;
    2557:           }
    2558:         }
    2559:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2560:           if (sieve.cone(*c_iter)->size() == 0) {
    2561:             renumbering[*c_iter] = max++;
    2562:           }
    2563:         }
    2564:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2565:           if (sieve.support(*b_iter)->size() == 0) {
    2566:             const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
    2567:             const typename Sieve::coneSequence::iterator coneEnd   = sieve.coneEnd(*b_iter);

    2569:             for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
    2570:               if (renumbering.find(*c_iter) == renumbering.end()) {
    2571:                 renumbering[*c_iter] = max++;
    2572:               }
    2573:             }
    2574:           }
    2575:         }
    2576:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2577:           if (renumbering.find(*b_iter) == renumbering.end()) {
    2578:             renumbering[*b_iter] = max++;
    2579:           }
    2580:         }
    2581:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2582:           if (renumbering.find(*c_iter) == renumbering.end()) {
    2583:             renumbering[*c_iter] = max++;
    2584:           }
    2585:         }
    2586:       } else {
    2587:         if (base->size()) {
    2588:           min = *base->begin();
    2589:           max = *base->begin();
    2590:         } else if (cap->size()) {
    2591:           min = *cap->begin();
    2592:           max = *cap->begin();
    2593:         }
    2594:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2595:           min = std::min(min, *b_iter);
    2596:           max = std::max(max, *b_iter);
    2597:         }
    2598:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2599:           min = std::min(min, *c_iter);
    2600:           max = std::max(max, *c_iter);
    2601:         }
    2602:         if (base->size() || cap->size()) {
    2603:           ++max;
    2604:         }
    2605:         for(typename ISieve::point_type p = min; p < max; ++p) {
    2606:           renumbering[p] = p;
    2607:         }
    2608:       }
    2609:       // Create the ISieve
    2610:       isieve.setChart(typename ISieve::chart_type(min, max));
    2611:       // Set cone and support sizes
    2612:       size_t maxSize = 0;

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

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

    2623:         isieve.setSupportSize(renumbering[*c_iter], support->size());
    2624:         maxSize = std::max(maxSize, support->size());
    2625:       }
    2626:       isieve.allocate();
    2627:       // Fill up cones and supports
    2628:       typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];

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

    2634:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
    2635:           points[i] = renumbering[*c_iter];
    2636:         }
    2637:         isieve.setCone(points, renumbering[*b_iter]);
    2638:       }
    2639:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2640:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
    2641:         int i = 0;

    2643:         for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
    2644:           points[i] = renumbering[*s_iter];
    2645:         }
    2646:         isieve.setSupport(renumbering[*c_iter], points);
    2647:       }
    2648:       delete [] points;
    2649:     }
    2650:     template<typename Sieve, typename Renumbering>
    2651:     static void convertSieve(Sieve& sieve, DM dm, Renumbering& renumbering, bool renumber = true) {
    2652:       // First construct a renumbering of the sieve points
    2653:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
    2654:       const Obj<typename Sieve::capSequence>&  cap  = sieve.cap();
    2655:       PetscInt                                 min  = 0;
    2656:       PetscInt                                 max  = 0;
    2657:       PetscErrorCode                           ierr;

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

    2662:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2663:           if (sieve.support(*b_iter)->size() == 0) {
    2664:             renumbering[*b_iter] = max++;
    2665:           }
    2666:         }
    2667:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2668:           if (sieve.cone(*c_iter)->size() == 0) {
    2669:             renumbering[*c_iter] = max++;
    2670:           }
    2671:         }
    2672: #if 0
    2673:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2674:           if (sieve.support(*b_iter)->size() == 0) {
    2675:             const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
    2676:             const typename Sieve::coneSequence::iterator coneEnd   = sieve.coneEnd(*b_iter);

    2678:             for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
    2679:               if (renumbering.find(*c_iter) == renumbering.end()) {
    2680:                 renumbering[*c_iter] = max++;
    2681:               }
    2682:             }
    2683:           }
    2684:         }
    2685: #else
    2686:         std::vector<typename Sieve::point_type> faces;
    2687:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2688:           if (sieve.support(*b_iter)->size() == 0) {
    2689:             const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
    2690:             const typename Sieve::coneSequence::iterator coneEnd   = sieve.coneEnd(*b_iter);

    2692:             for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
    2693:               if (renumbering.find(*c_iter) == renumbering.end()) {
    2694:                 faces.push_back(*c_iter);
    2695:               }
    2696:             }
    2697:           }
    2698:         }
    2699:         std::sort(faces.begin(), faces.end());
    2700:         typename std::vector<typename Sieve::point_type>::const_iterator fEnd = std::unique(faces.begin(), faces.end());
    2701:         for(typename std::vector<typename Sieve::point_type>::const_iterator c_iter = faces.begin(); c_iter != fEnd; ++c_iter) {
    2702:           renumbering[*c_iter] = max++;
    2703:         }
    2704:         faces.clear();
    2705: #endif
    2706:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2707:           if (renumbering.find(*b_iter) == renumbering.end()) {
    2708:             renumbering[*b_iter] = max++;
    2709:           }
    2710:         }
    2711:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2712:           if (renumbering.find(*c_iter) == renumbering.end()) {
    2713:             renumbering[*c_iter] = max++;
    2714:           }
    2715:         }
    2716:       } else {
    2717:         if (base->size()) {
    2718:           min = *base->begin();
    2719:           max = *base->begin();
    2720:         } else if (cap->size()) {
    2721:           min = *cap->begin();
    2722:           max = *cap->begin();
    2723:         }
    2724:         for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
    2725:           min = std::min(min, *b_iter);
    2726:           max = std::max(max, *b_iter);
    2727:         }
    2728:         for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2729:           min = std::min(min, *c_iter);
    2730:           max = std::max(max, *c_iter);
    2731:         }
    2732:         if (base->size() || cap->size()) {
    2733:           ++max;
    2734:         }
    2735:         for(PetscInt p = min; p < max; ++p) {
    2736:           renumbering[p] = p;
    2737:         }
    2738:       }
    2739:       // Create the ISieve
    2740:       DMPlexSetChart(dm, min, max);CHKERRXX(ierr);
    2741:       // Set cone and support sizes
    2742:       size_t maxSize = 0;

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

    2747:         DMPlexSetConeSize(dm, renumbering[*b_iter], cone->size());CHKERRXX(ierr);
    2748:         maxSize = std::max(maxSize, cone->size());
    2749:       }
    2750:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2751:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);

    2753:         DMPlexSetSupportSize(dm, renumbering[*c_iter], support->size());CHKERRXX(ierr);
    2754:         maxSize = std::max(maxSize, support->size());
    2755:       }
    2756:       DMSetUp(dm);CHKERRXX(ierr);
    2757:       // Fill up cones and supports
    2758:       typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];

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

    2764:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
    2765:           points[i] = renumbering[*c_iter];
    2766:         }
    2767:         DMPlexSetCone(dm, renumbering[*b_iter], points);CHKERRXX(ierr);
    2768:       }
    2769:       for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
    2770:         const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
    2771:         int i = 0;

    2773:         for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
    2774:           points[i] = renumbering[*s_iter];
    2775:         }
    2776:         DMPlexSetSupport(dm, renumbering[*c_iter], points);CHKERRXX(ierr);
    2777:       }
    2778:       delete [] points;
    2779:     }
    2780:     template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
    2781:     static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
    2782:       if (isieve.getMaxConeSize() < 0) return;
    2783:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
    2784:       int *orientations = new int[isieve.getMaxConeSize()];

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

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

    2793:           orientations[i] = orientation->restrictPoint(arrow)[0];
    2794:         }
    2795:         isieve.setConeOrientation(orientations, renumbering[*b_iter]);
    2796:       }
    2797:       delete [] orientations;
    2798:     }
    2799:     template<typename Sieve, typename Renumbering, typename ArrowSection>
    2800:     static void convertOrientation(Sieve& sieve, DM dm, Renumbering& renumbering, ArrowSection *orientation) {
    2801:       PetscInt       maxConeSize;

    2804:       DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRXX(ierr);
    2805:       if (maxConeSize < 0) return;
    2806:       const Obj<typename Sieve::baseSequence>& base = sieve.base();
    2807:       int *orientations = new int[maxConeSize];

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

    2813:         for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
    2814:           typename ArrowSection::point_type arrow(*c_iter, *b_iter);
    2815:           const int o = orientation->restrictPoint(arrow)[0];

    2817:           orientations[i] = o == 1 ? 0 : o;
    2818:         }
    2819:         DMPlexSetConeOrientation(dm, renumbering[*b_iter], orientations);
    2820:       }
    2821:       delete [] orientations;
    2822:     }
    2823:     template<typename Section, typename ISection, typename Renumbering>
    2824:     static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
    2825:       const typename Section::chart_type& chart = coordinates.getChart();
    2826:       typename ISection::point_type       min   = *chart.begin();
    2827:       typename ISection::point_type       max   = *chart.begin();

    2829:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2830:         min = std::min(min, renumbering[*p_iter]);
    2831:         max = std::max(max, renumbering[*p_iter]);
    2832:       }
    2833:       icoordinates.setChart(typename ISection::chart_type(min, max+1));
    2834:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2835:         icoordinates.setFiberDimension(renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));
    2836:       }
    2837:       icoordinates.allocatePoint();
    2838:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2839:         icoordinates.updatePoint(renumbering[*p_iter], coordinates.restrictPoint(*p_iter));
    2840:       }
    2841:     }
    2842:     template<typename Section, typename Renumbering>
    2843:     static void convertCoordinates(Section& coordinates, PetscSection coordSection, Vec coords, Renumbering& renumbering) {
    2844:       const typename Section::chart_type& chart = coordinates.getChart();
    2845:       PetscInt                            min   = *chart.begin();
    2846:       PetscInt                            max   = *chart.begin();
    2847:       PetscScalar                        *a;
    2848:       PetscInt                            n;
    2849:       PetscErrorCode                      ierr;

    2851:       PetscSectionSetNumFields(coordSection, 1);CHKERRXX(ierr);
    2852:       if (!chart.size()) {
    2853:         PetscSectionSetFieldComponents(coordSection, 0, 1);CHKERRXX(ierr);
    2854:         return;
    2855:       }
    2856:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2857:         min = std::min(min, renumbering[*p_iter]);
    2858:         max = std::max(max, renumbering[*p_iter]);
    2859:       }
    2860:       PetscSectionSetFieldComponents(coordSection, 0, coordinates.getFiberDimension(*chart.begin()));CHKERRXX(ierr);
    2861:       PetscSectionSetChart(coordSection, min, max+1);CHKERRXX(ierr);
    2862:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2863:         PetscSectionSetDof(coordSection, renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
    2864:         PetscSectionSetFieldDof(coordSection, renumbering[*p_iter], 0, coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
    2865:       }
    2866:       PetscSectionSetUp(coordSection);CHKERRXX(ierr);
    2867:       PetscSectionGetStorageSize(coordSection, &n);CHKERRXX(ierr);
    2868:       VecSetSizes(coords, n, PETSC_DETERMINE);CHKERRXX(ierr);
    2869:       VecSetFromOptions(coords);CHKERRXX(ierr);
    2870:       VecGetArray(coords, &a);CHKERRXX(ierr);
    2871:       for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
    2872:         const typename Section::value_type *values = coordinates.restrictPoint(*p_iter);
    2873:         PetscInt dof, off;

    2875:         PetscSectionGetDof(coordSection, renumbering[*p_iter], &dof);CHKERRXX(ierr);
    2876:         PetscSectionGetOffset(coordSection, renumbering[*p_iter], &off);CHKERRXX(ierr);
    2877:         for(int d = 0; d < dof; ++d) {
    2878:           a[off+d] = values[d];
    2879:         }
    2880:       }
    2881:       VecRestoreArray(coords, &a);CHKERRXX(ierr);
    2882:     }
    2883:     template<typename Label, typename Renumbering>
    2884:     static void convertLabel(const Obj<Label>& newLabel, const Obj<Label>& oldLabel, Renumbering& renumbering) {
    2885:       for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
    2886:         if (oldLabel->getConeSize(p->first)) {
    2887:           newLabel->setCone(*oldLabel->cone(p->first)->begin(), p->second);
    2888:         }
    2889:       }
    2890:     }
    2891:     template<typename Label, typename Renumbering>
    2892:     static void convertLabel(DM dm, const char name[], const Obj<Label>& label, Renumbering& renumbering) {

    2895:       for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
    2896:         if (label->getConeSize(p->first)) {
    2897:           DMPlexSetLabelValue(dm, name, p->second, *label->cone(p->first)->begin());CHKERRXX(ierr);
    2898:         }
    2899:       }
    2900:     }
    2901:     template<typename Mesh, typename IMesh, typename Renumbering>
    2902:     static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
    2903:       convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
    2904:       imesh.stratify();
    2905:       convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
    2906:       convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
    2907:       if (mesh.hasRealSection("normals")) {
    2908:         convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
    2909:       }
    2910:       const typename Mesh::labels_type& labels = mesh.getLabels();
    2911:       std::string heightName("height");
    2912:       std::string depthName("depth");

    2914:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
    2915: #ifdef IMESH_NEW_LABELS
    2916:         if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
    2917:           convertLabel(imesh, l_iter->first, l_iter->second);
    2918:         }
    2919: #else
    2920:         if (renumber) {
    2921:           convertLabel(imesh.createLabel(l_iter->first), l_iter->second, renumbering);
    2922:         } else {
    2923:           imesh.setLabel(l_iter->first, l_iter->second);
    2924:         }
    2925: #endif
    2926:       }
    2927:     }
    2928:     template<typename Mesh, typename Renumbering>
    2929:     static void convertMesh(Mesh& mesh, DM *dm, Renumbering& renumbering, bool renumber = true) {
    2930:       PetscSection   coordSection;
    2931:       Vec            coordinates;

    2934:       DMCreate(mesh.comm(), dm);CHKERRXX(ierr);
    2935:       DMSetType(*dm, DMPLEX);CHKERRXX(ierr);
    2936:       DMPlexSetDimension(*dm, mesh.getDimension());CHKERRXX(ierr);
    2937:       convertSieve(*mesh.getSieve(), *dm, renumbering, renumber);
    2938:       DMPlexStratify(*dm);CHKERRXX(ierr);
    2939:       convertOrientation(*mesh.getSieve(), *dm, renumbering, mesh.getArrowSection("orientation").ptr());
    2940:       DMPlexGetCoordinateSection(*dm, &coordSection);CHKERRXX(ierr);
    2941:       VecCreate(mesh.comm(), &coordinates);CHKERRXX(ierr);
    2942:       convertCoordinates(*mesh.getRealSection("coordinates"), coordSection, coordinates, renumbering);
    2943:       DMSetCoordinatesLocal(*dm, coordinates);CHKERRXX(ierr);
    2944:       VecDestroy(&coordinates);CHKERRXX(ierr);
    2945:       const typename Mesh::labels_type& labels = mesh.getLabels();

    2947:       for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
    2948:         convertLabel(*dm, l_iter->first.c_str(), l_iter->second, renumbering);
    2949:       }
    2950:     }
    2951:   };

    2953:   class ISieveSerializer {
    2954:   public:
    2955:     template<typename ISieve>
    2956:     static void writeSieve(const std::string& filename, ISieve& sieve) {
    2957:       std::ofstream fs;

    2959:       if (sieve.commRank() == 0) {
    2960:         fs.open(filename.c_str());
    2961:       }
    2962:       writeSieve(fs, sieve);
    2963:       if (sieve.commRank() == 0) {
    2964:         fs.close();
    2965:       }
    2966:     }
    2967:     template<typename ISieve>
    2968:     static void writeSieve(std::ofstream& fs, ISieve& sieve) {
    2969:       typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
    2970:       const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
    2971:       typename ISieve::point_type             min   = chart->min();
    2972:       typename ISieve::point_type             max   = chart->max();
    2973:       PetscInt                               *mins  = new PetscInt[sieve.commSize()];
    2974:       PetscInt                               *maxs  = new PetscInt[sieve.commSize()];

    2976:       // Write sizes
    2977:       if (sieve.commRank() == 0) {
    2978:         // Write local
    2979:         fs << min <<" "<< max << std::endl;
    2980:         for(typename ISieve::point_type p = min; p < max; ++p) {
    2981:           fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
    2982:         }
    2983:         // Receive and write remote
    2984:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
    2985:           PetscInt       minp, maxp;
    2986:           PetscInt      *sizes;
    2987:           MPI_Status     status;

    2990:           MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
    2991:           MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
    2992:           PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
    2993:           MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
    2994:           fs << minp <<" "<< maxp << std::endl;
    2995:           for(PetscInt s = 0; s < maxp - minp; ++s) {
    2996:             fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
    2997:           }
    2998:           PetscFree(sizes);CHKERRXX(ierr);
    2999:           mins[pr] = minp;
    3000:           maxs[pr] = maxp;
    3001:         }
    3002:       } else {
    3003:         // Send remote
    3004:         //PetscInt       min = chart->min();
    3005:         //PetscInt       max = chart->max();
    3006:         PetscInt       s   = 0;
    3007:         PetscInt      *sizes;

    3010:         MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
    3011:         MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
    3012:         PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
    3013:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
    3014:           sizes[s*2+0] = sieve.getConeSize(p);
    3015:           sizes[s*2+1] = sieve.getSupportSize(p);
    3016:         }
    3017:         MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
    3018:         PetscFree(sizes);CHKERRXX(ierr);
    3019:       }
    3020:       // Write covering
    3021:       if (sieve.commRank() == 0) {
    3022:         // Write local
    3023:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);

    3025:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3026:           sieve.cone(p, pV);
    3027:           const typename Visitor::point_type *cone  = pV.getPoints();
    3028:           const int                           cSize = pV.getSize();

    3030:           if (cSize > 0) {
    3031:             for(int c = 0; c < cSize; ++c) {
    3032:               if (c) {fs << " ";}
    3033:               fs << cone[c];
    3034:             }
    3035:             fs << std::endl;
    3036:           }
    3037:           pV.clear();

    3039:           sieve.orientedCone(p, pV);
    3040:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
    3041:           const int                                    oSize = pV.getOrientedSize();

    3043:           if (oSize > 0) {
    3044:             for(int c = 0; c < oSize; ++c) {
    3045:               if (c) {fs << " ";}
    3046:               fs << oCone[c].second;
    3047:             }
    3048:             fs << std::endl;
    3049:           }
    3050:           pV.clear();

    3052:           sieve.support(p, pV);
    3053:           const typename Visitor::point_type *support = pV.getPoints();
    3054:           const int                           sSize   = pV.getSize();

    3056:           if (sSize > 0) {
    3057:             for(int s = 0; s < sSize; ++s) {
    3058:               if (s) {fs << " ";}
    3059:               fs << support[s];
    3060:             }
    3061:             fs << std::endl;
    3062:           }
    3063:           pV.clear();
    3064:         }
    3065:         // Receive and write remote
    3066:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
    3067:           PetscInt       size;
    3068:           PetscInt      *data;
    3069:           PetscInt       off = 0;
    3070:           MPI_Status     status;

    3073:           MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3074:           PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
    3075:           MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3076:           for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
    3077:             PetscInt cSize = data[off++];

    3079:             fs << cSize << std::endl;
    3080:             if (cSize > 0) {
    3081:               for(int c = 0; c < cSize; ++c) {
    3082:                 if (c) {fs << " ";}
    3083:                 fs << data[off++];
    3084:               }
    3085:               fs << std::endl;
    3086:             }
    3087:             PetscInt oSize = data[off++];

    3089:             if (oSize > 0) {
    3090:               for(int c = 0; c < oSize; ++c) {
    3091:                 if (c) {fs << " ";}
    3092:                 fs << data[off++];
    3093:               }
    3094:               fs << std::endl;
    3095:             }
    3096:             PetscInt sSize = data[off++];

    3098:             fs << sSize << std::endl;
    3099:             if (sSize > 0) {
    3100:               for(int s = 0; s < sSize; ++s) {
    3101:                 if (s) {fs << " ";}
    3102:                 fs << data[off++];
    3103:               }
    3104:               fs << std::endl;
    3105:             }
    3106:           }
    3107:           assert(off == size);
    3108:           PetscFree(data);CHKERRXX(ierr);
    3109:         }
    3110:       } else {
    3111:         // Send remote
    3112:         Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
    3113:         PetscInt totalConeSize    = 0;
    3114:         PetscInt totalSupportSize = 0;

    3116:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3117:           totalConeSize    += sieve.getConeSize(p);
    3118:           totalSupportSize += sieve.getSupportSize(p);
    3119:         }
    3120:         PetscInt       size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
    3121:         PetscInt       off  = 0;
    3122:         PetscInt      *data;

    3125:         MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
    3126:         // There is no nice way to make a generic MPI type here. Really sucky
    3127:         PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
    3128:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3129:           sieve.cone(p, pV);
    3130:           const typename Visitor::point_type *cone  = pV.getPoints();
    3131:           const int                           cSize = pV.getSize();

    3133:           data[off++] = cSize;
    3134:           for(int c = 0; c < cSize; ++c) {
    3135:             data[off++] = cone[c];
    3136:           }
    3137:           pV.clear();

    3139:           sieve.orientedCone(p, pV);
    3140:           const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
    3141:           const int                                    oSize = pV.getOrientedSize();

    3143:           data[off++] = oSize;
    3144:           for(int c = 0; c < oSize; ++c) {
    3145:             data[off++] = oCone[c].second;
    3146:           }
    3147:           pV.clear();

    3149:           sieve.support(p, pV);
    3150:           const typename Visitor::point_type *support = pV.getPoints();
    3151:           const int                           sSize   = pV.getSize();

    3153:           data[off++] = sSize;
    3154:           for(int s = 0; s < sSize; ++s) {
    3155:             data[off++] = support[s];
    3156:           }
    3157:           pV.clear();
    3158:         }
    3159:         MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
    3160:         PetscFree(data);CHKERRXX(ierr);
    3161:       }
    3162:       delete [] mins;
    3163:       delete [] maxs;
    3164:       // Output renumbering
    3165:     }
    3166:     template<typename ISieve>
    3167:     static void loadSieve(const std::string& filename, ISieve& sieve) {
    3168:       std::ifstream fs;

    3170:       if (sieve.commRank() == 0) {
    3171:         fs.open(filename.c_str());
    3172:       }
    3173:       loadSieve(fs, sieve);
    3174:       if (sieve.commRank() == 0) {
    3175:         fs.close();
    3176:       }
    3177:     }
    3178:     template<typename ISieve>
    3179:     static void loadSieve(std::ifstream& fs, ISieve& sieve) {
    3180:       typename ISieve::point_type min, max;
    3181:       PetscInt                   *mins = new PetscInt[sieve.commSize()];
    3182:       PetscInt                   *maxs = new PetscInt[sieve.commSize()];
    3183:       PetscInt                   *totalConeSizes    = new PetscInt[sieve.commSize()];
    3184:       PetscInt                   *totalSupportSizes = new PetscInt[sieve.commSize()];

    3186:       // Load sizes
    3187:       if (sieve.commRank() == 0) {
    3188:         // Load local
    3189:         fs >> min;
    3190:         fs >> max;
    3191:         sieve.setChart(typename ISieve::chart_type(min, max));
    3192:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3193:           typename ISieve::index_type coneSize, supportSize;

    3195:           fs >> coneSize;
    3196:           fs >> supportSize;
    3197:           sieve.setConeSize(p, coneSize);
    3198:           sieve.setSupportSize(p, supportSize);
    3199:         }
    3200:         // Load and send remote
    3201:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
    3202:           PetscInt       minp, maxp;
    3203:           PetscInt      *sizes;

    3206:           fs >> minp;
    3207:           fs >> maxp;
    3208:           MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
    3209:           MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
    3210:           PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
    3211:           totalConeSizes[pr]    = 0;
    3212:           totalSupportSizes[pr] = 0;
    3213:           for(PetscInt s = 0; s < maxp - minp; ++s) {
    3214:             fs >> sizes[s*2+0];
    3215:             fs >> sizes[s*2+1];
    3216:             totalConeSizes[pr]    += sizes[s*2+0];
    3217:             totalSupportSizes[pr] += sizes[s*2+1];
    3218:           }
    3219:           MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
    3220:           PetscFree(sizes);CHKERRXX(ierr);
    3221:           mins[pr] = minp;
    3222:           maxs[pr] = maxp;
    3223:         }
    3224:       } else {
    3225:         // Load remote
    3226:         PetscInt       s   = 0;
    3227:         PetscInt      *sizes;
    3228:         MPI_Status     status;

    3231:         MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3232:         MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3233:         sieve.setChart(typename ISieve::chart_type(min, max));
    3234:         PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
    3235:         MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3236:         for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
    3237:           sieve.setConeSize(p, sizes[s*2+0]);
    3238:           sieve.setSupportSize(p, sizes[s*2+1]);
    3239:         }
    3240:         PetscFree(sizes);CHKERRXX(ierr);
    3241:       }
    3242:       sieve.allocate();
    3243:       // Load covering
    3244:       if (sieve.commRank() == 0) {
    3245:         // Load local
    3246:         typename ISieve::index_type  maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
    3247:         typename ISieve::point_type *points  = new typename ISieve::point_type[maxSize];

    3249:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3250:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
    3251:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);

    3253:           if (coneSize > 0) {
    3254:             for(int c = 0; c < coneSize; ++c) {
    3255:               fs >> points[c];
    3256:             }
    3257:             sieve.setCone(points, p);
    3258:             if (sieve.orientedCones()) {
    3259:               for(int c = 0; c < coneSize; ++c) {
    3260:                 fs >> points[c];
    3261:               }
    3262:               sieve.setConeOrientation(points, p);
    3263:             }
    3264:           }
    3265:           if (supportSize > 0) {
    3266:             for(int s = 0; s < supportSize; ++s) {
    3267:               fs >> points[s];
    3268:             }
    3269:             sieve.setSupport(p, points);
    3270:           }
    3271:         }
    3272:         delete [] points;
    3273:         // Load and send remote
    3274:         for(int pr = 1; pr < sieve.commSize(); ++pr) {
    3275:           PetscInt       size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
    3276:           PetscInt       off  = 0;
    3277:           PetscInt      *data;

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

    3286:             fs >> coneSize;
    3287:             data[off++] = coneSize;
    3288:             if (coneSize > 0) {
    3289:               for(int c = 0; c < coneSize; ++c) {
    3290:                 fs >> data[off++];
    3291:               }
    3292:               for(int c = 0; c < coneSize; ++c) {
    3293:                 fs >> data[off++];
    3294:               }
    3295:             }
    3296:             fs >> supportSize;
    3297:             data[off++] = supportSize;
    3298:             if (supportSize > 0) {
    3299:               for(int s = 0; s < supportSize; ++s) {
    3300:                 fs >> data[off++];
    3301:               }
    3302:             }
    3303:           }
    3304:           assert(off == size);
    3305:           MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
    3306:           PetscFree(data);CHKERRXX(ierr);
    3307:         }
    3308:         delete [] mins;
    3309:         delete [] maxs;
    3310:       } else {
    3311:         // Load remote
    3312:         PetscInt       size;
    3313:         PetscInt      *data;
    3314:         PetscInt       off = 0;
    3315:         MPI_Status     status;

    3318:         MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3319:         PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
    3320:         MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
    3321:         for(typename ISieve::point_type p = min; p < max; ++p) {
    3322:           typename ISieve::index_type coneSize    = sieve.getConeSize(p);
    3323:           typename ISieve::index_type supportSize = sieve.getSupportSize(p);
    3324:           PetscInt cs = data[off++];

    3326:           assert(cs == coneSize);
    3327:           if (coneSize > 0) {
    3328:             sieve.setCone(&data[off], p);
    3329:             off += coneSize;
    3330:             if (sieve.orientedCones()) {
    3331:               sieve.setConeOrientation(&data[off], p);
    3332:               off += coneSize;
    3333:             }
    3334:           }
    3335:           PetscInt ss = data[off++];

    3337:           assert(ss == supportSize);
    3338:           if (supportSize > 0) {
    3339:             sieve.setSupport(p, &data[off]);
    3340:             off += supportSize;
    3341:           }
    3342:         }
    3343:         assert(off == size);
    3344:       }
    3345:       // Load renumbering
    3346:     }
    3347:   };
    3348: }

    3350: #endif