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

      1: #ifndef included_ALE_Sections_hh
      2: #define included_ALE_Sections_hh

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

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

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

    106:       if (size > maxSize) {
    107:         const value_type dummy(0);

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

    136:       for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    137:         array[c++] = *c_iter;
    138:       }
    139:       return array;
    140:     };
    141:   };

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

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

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

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

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

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

    398:   template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
    399:   class LabelBaseSectionV : public ALE::ParallelObject {
    400:   public:
    401:     typedef Sieve_                                    sieve_type;
    402:     typedef Label_                                    label_type;
    403:     typedef Alloc_                                    alloc_type;
    404:     typedef int                                       value_type;
    405:     typedef typename sieve_type::target_type          point_type;
    406:     //typedef typename sieve_type::traits::baseSequence chart_type;
    407:     typedef int                                       chart_type;
    408:   protected:
    409:     Obj<sieve_type> _sieve;
    410:     Obj<label_type> _label;
    411:     //chart_type      _chart;
    412:     int             _sizes[2];
    413:   public:
    414:     //LabelBaseSectionV(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
    415:     LabelBaseSectionV(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label) {_sizes[0] = 1; _sizes[1] = 0;};
    416:     ~LabelBaseSectionV() {};
    417:   public: // Verifiers
    418:     bool hasPoint(const point_type& point) const {
    419:       return this->_label->cone(point)->size() ? true : false;
    420:     };
    421:   public:
    422:     //const chart_type& getChart() const {
    423:     //  return this->_chart;
    424:     //};
    425:     int getFiberDimension(const point_type& p) const {
    426:       return this->hasPoint(p) ? 1 : 0;
    427:     };
    428:     const value_type *restrictSpace() const {
    429:       return this->_sizes;
    430:     };
    431:     const value_type *restrictPoint(const point_type& p) const {
    432:       if (this->hasPoint(p)) return this->_sizes;
    433:       return &this->_sizes[1];
    434:     };
    435:   };

    437:   namespace New {
    438:     // This section takes an existing section, and reports instead the fiber dimensions as values
    439:     template<typename Section_>
    440:     class SizeSection : public ALE::ParallelObject {
    441:     public:
    442:       typedef Section_                          section_type;
    443:       typedef typename section_type::point_type point_type;
    444:       typedef int                               value_type;
    445:     protected:
    446:       Obj<section_type> _section;
    447:       value_type        _size;
    448:     public:
    449:       SizeSection(const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, section->debug()), _section(section) {};
    450:       virtual ~SizeSection() {};
    451:     public:
    452:       bool hasPoint(const point_type& point) {
    453:         return this->_section->hasPoint(point);
    454:       };
    455:       const value_type *restrictPoint(const point_type& p) {
    456:         this->_size = this->_section->getFiberDimension(p);
    457:         return &this->_size;
    458:       };
    459:     public:
    460:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    461:         this->_section->view(name, comm);
    462:       };
    463:     };

    465:     // This section reports as values the size of the partition associated with the partition point
    466:     template<typename Bundle_, typename Marker_>
    467:     class PartitionSizeSection : public ALE::ParallelObject {
    468:     public:
    469:       typedef Bundle_                          bundle_type;
    470:       typedef typename bundle_type::sieve_type sieve_type;
    471:       typedef typename bundle_type::point_type point_type;
    472:       typedef Marker_                          marker_type;
    473:       typedef int                              value_type;
    474:       typedef std::map<marker_type, int>       sizes_type;
    475:     protected:
    476:       sizes_type _sizes;
    477:       int        _height;
    478:       void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
    479:         const Obj<typename bundle_type::label_sequence>& cells      = bundle->heightStratum(this->_height);
    480:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
    481:         std::map<marker_type, std::set<point_type> >     points;

    483:         if (numElements != (int) cells->size()) {
    484:           throw ALE::Exception("Partition size does not match the number of elements");
    485:         }
    486:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
    487:           typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
    488:           const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
    489:           const int idx = cNumbering->getIndex(*e_iter);

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

    495:             points[partition[idx]].insert(star->begin(), star->end());
    496:           }
    497:         }
    498:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    499:           this->_sizes[p_iter->first] = p_iter->second.size();
    500:         }
    501:       };
    502:     public:
    503:       PartitionSizeSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
    504:         this->_init(bundle, numElements, partition);
    505:       };
    506:       virtual ~PartitionSizeSection() {};
    507:     public:
    508:       bool hasPoint(const point_type& point) {return true;};
    509:       const value_type *restrictPoint(const point_type& p) {
    510:         return &this->_sizes[p];
    511:       };
    512:     public:
    513:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    514:         ostringstream txt;
    515:         int rank;

    517:         if (comm == MPI_COMM_NULL) {
    518:           comm = this->comm();
    519:           rank = this->commRank();
    520:         } else {
    521:           MPI_Comm_rank(comm, &rank);
    522:         }
    523:         if (name == "") {
    524:           if(rank == 0) {
    525:             txt << "viewing a PartitionSizeSection" << std::endl;
    526:           }
    527:         } else {
    528:           if(rank == 0) {
    529:             txt << "viewing PartitionSizeSection '" << name << "'" << std::endl;
    530:           }
    531:         }
    532:         for(typename sizes_type::const_iterator s_iter = this->_sizes.begin(); s_iter != this->_sizes.end(); ++s_iter) {
    533:           const marker_type& partition = s_iter->first;
    534:           const value_type   size      = s_iter->second;

    536:           txt << "[" << this->commRank() << "]: Partition " << partition << " size " << size << std::endl;
    537:         }
    538:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    539:         PetscSynchronizedFlush(comm);
    540:       };
    541:     };

    543:     template<typename Point_>
    544:     class PartitionDomain {
    545:     public:
    546:       typedef Point_ point_type;
    547:     public:
    548:       PartitionDomain() {};
    549:       ~PartitionDomain() {};
    550:     public:
    551:       int count(const point_type& point) const {return 1;};
    552:     };

    554:     // This section returns the points in each partition
    555:     template<typename Bundle_, typename Marker_>
    556:     class PartitionSection : public ALE::ParallelObject {
    557:     public:
    558:       typedef Bundle_                            bundle_type;
    559:       typedef typename bundle_type::sieve_type   sieve_type;
    560:       typedef typename bundle_type::point_type   point_type;
    561:       typedef Marker_                            marker_type;
    562:       typedef int                                value_type;
    563:       typedef std::map<marker_type, point_type*> points_type;
    564:       typedef PartitionDomain<point_type>        chart_type;
    565:     protected:
    566:       points_type _points;
    567:       chart_type  _domain;
    568:       int         _height;
    569:       void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
    570:         // Should check for patch 0
    571:         const Obj<typename bundle_type::label_sequence>& cells      = bundle->heightStratum(this->_height);
    572:         const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
    573:         std::map<marker_type, std::set<point_type> >     points;
    574:         std::map<marker_type, int>                       offsets;

    576:         if (numElements != (int) cells->size()) {
    577:           throw ALE::Exception("Partition size does not match the number of elements");
    578:         }
    579:         for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
    580:           typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
    581:           const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
    582:           const int idx = cNumbering->getIndex(*e_iter);

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

    588:             points[partition[idx]].insert(star->begin(), star->end());
    589:           }
    590:         }
    591:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    592:           this->_points[p_iter->first] = new point_type[p_iter->second.size()];
    593:           offsets[p_iter->first] = 0;
    594:           for(typename std::set<point_type>::const_iterator s_iter = p_iter->second.begin(); s_iter != p_iter->second.end(); ++s_iter) {
    595:             this->_points[p_iter->first][offsets[p_iter->first]++] = *s_iter;
    596:           }
    597:         }
    598:         for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
    599:           if (offsets[p_iter->first] != (int) p_iter->second.size()) {
    600:             ostringstream txt;
    601:             txt << "Invalid offset for partition " << p_iter->first << ": " << offsets[p_iter->first] << " should be " << p_iter->second.size();
    602:             throw ALE::Exception(txt.str().c_str());
    603:           }
    604:         }
    605:       };
    606:     public:
    607:       PartitionSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
    608:         this->_init(bundle, numElements, partition);
    609:       };
    610:       virtual ~PartitionSection() {
    611:         for(typename points_type::iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
    612:           delete [] p_iter->second;
    613:         }
    614:       };
    615:     public:
    616:       const chart_type& getChart() {return this->_domain;};
    617:       bool hasPoint(const point_type& point) {return true;};
    618:       const value_type *restrictPoint(const point_type& p) {
    619:         return this->_points[p];
    620:       };
    621:     public:
    622:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    623:         ostringstream txt;
    624:         int rank;

    626:         if (comm == MPI_COMM_NULL) {
    627:           comm = this->comm();
    628:           rank = this->commRank();
    629:         } else {
    630:           MPI_Comm_rank(comm, &rank);
    631:         }
    632:         if (name == "") {
    633:           if(rank == 0) {
    634:             txt << "viewing a PartitionSection" << std::endl;
    635:           }
    636:         } else {
    637:           if(rank == 0) {
    638:             txt << "viewing PartitionSection '" << name << "'" << std::endl;
    639:           }
    640:         }
    641:         for(typename points_type::const_iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
    642:           const marker_type& partition  = p_iter->first;
    643:           //const point_type *points = p_iter->second;

    645:           txt << "[" << this->commRank() << "]: Partition " << partition << std::endl;
    646:         }
    647:         if (this->_points.size() == 0) {
    648:           txt << "[" << this->commRank() << "]: empty" << std::endl;
    649:         }
    650:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    651:         PetscSynchronizedFlush(comm);
    652:       };
    653:     };

    655:     template<typename Bundle_, typename Sieve_>
    656:     class ConeSizeSection : public ALE::ParallelObject {
    657:     public:
    658:       typedef ConeSizeSection<Bundle_, Sieve_> section_type;
    659:       typedef int                              patch_type;
    660:       typedef Bundle_                          bundle_type;
    661:       typedef Sieve_                           sieve_type;
    662:       typedef typename bundle_type::point_type point_type;
    663:       typedef int                              value_type;
    664:     protected:
    665:       Obj<bundle_type> _bundle;
    666:       Obj<sieve_type>  _sieve;
    667:       value_type       _size;
    668:       int              _minHeight;
    669:       Obj<section_type> _section;
    670:     public:
    671:       ConeSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumHeight = 0) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _bundle(bundle), _sieve(sieve), _minHeight(minimumHeight) {
    672:         this->_section = this;
    673:         this->_section.addRef();
    674:       };
    675:       virtual ~ConeSizeSection() {};
    676:     public: // Verifiers
    677:       bool hasPoint(const point_type& point) {return true;};
    678:     public: // Restriction
    679:       const value_type *restrictPoint(const point_type& p) {
    680:         if ((this->_minHeight == 0) || (this->_bundle->height(p) >= this->_minHeight)) {
    681:           this->_size = this->_sieve->cone(p)->size();
    682:         } else {
    683:           this->_size = 0;
    684:         }
    685:         return &this->_size;
    686:       };
    687:     public: // Adapter
    688:       const Obj<section_type>& getSection(const patch_type& patch) {
    689:         return this->_section;
    690:       };
    691:     public:
    692:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    693:         ostringstream txt;
    694:         int rank;

    696:         if (comm == MPI_COMM_NULL) {
    697:           comm = this->comm();
    698:           rank = this->commRank();
    699:         } else {
    700:           MPI_Comm_rank(comm, &rank);
    701:         }
    702:         if (name == "") {
    703:           if(rank == 0) {
    704:             txt << "viewing a ConeSizeSection" << std::endl;
    705:           }
    706:         } else {
    707:           if(rank == 0) {
    708:             txt << "viewing ConeSizeSection '" << name << "'" << std::endl;
    709:           }
    710:         }
    711:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    712:         PetscSynchronizedFlush(comm);
    713:       };
    714:     };

    716:     template<typename Sieve_>
    717:     class ConeSection : public ALE::ParallelObject {
    718:     public:
    719:       typedef Sieve_                           sieve_type;
    720:       typedef typename sieve_type::target_type point_type;
    721:       typedef typename sieve_type::source_type value_type;
    722:       typedef PartitionDomain<sieve_type>      chart_type;
    723:     protected:
    724:       Obj<sieve_type> _sieve;
    725:       int             _coneSize;
    726:       value_type     *_cone;
    727:       chart_type      _domain;
    728:       void ensureCone(const int size) {
    729:         if (size > this->_coneSize) {
    730:           if (this->_cone) delete [] this->_cone;
    731:           this->_coneSize = size;
    732:           this->_cone     = new value_type[this->_coneSize];
    733:         }
    734:       };
    735:     public:
    736:       ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _coneSize(-1), _cone(NULL) {};
    737:       virtual ~ConeSection() {if (this->_cone) delete [] this->_cone;};
    738:     public:
    739:       const chart_type& getChart() {return this->_domain;};
    740:       bool hasPoint(const point_type& point) {return true;};
    741:       const value_type *restrictPoint(const point_type& p) {
    742:         const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
    743:         int c = 0;

    745:         this->ensureCone(cone->size());
    746:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    747:           this->_cone[c++] = *c_iter;
    748:         }
    749:         return this->_cone;
    750:       };
    751:     public:
    752:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    753:         ostringstream txt;
    754:         int rank;

    756:         if (comm == MPI_COMM_NULL) {
    757:           comm = this->comm();
    758:           rank = this->commRank();
    759:         } else {
    760:           MPI_Comm_rank(comm, &rank);
    761:         }
    762:         if (name == "") {
    763:           if(rank == 0) {
    764:             txt << "viewing a ConeSection" << std::endl;
    765:           }
    766:         } else {
    767:           if(rank == 0) {
    768:             txt << "viewing ConeSection '" << name << "'" << std::endl;
    769:           }
    770:         }
    771:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    772:         PetscSynchronizedFlush(comm);
    773:       };
    774:     };

    776:     template<typename Bundle_, typename Sieve_>
    777:     class SupportSizeSection : public ALE::ParallelObject {
    778:     public:
    779:       typedef Bundle_                          bundle_type;
    780:       typedef Sieve_                           sieve_type;
    781:       typedef typename sieve_type::source_type point_type;
    782:       typedef typename sieve_type::target_type value_type;
    783:     protected:
    784:       Obj<bundle_type> _bundle;
    785:       Obj<sieve_type>  _sieve;
    786:       value_type       _size;
    787:       int              _minDepth;
    788:     public:
    789:       SupportSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumDepth = 0) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _bundle(bundle), _sieve(sieve), _minDepth(minimumDepth) {};
    790:       virtual ~SupportSizeSection() {};
    791:     public:
    792:       bool hasPoint(const point_type& point) {return true;};
    793:       const value_type *restrictPoint(const point_type& p) {
    794:         if ((this->_minDepth == 0) || (this->_bundle->depth(p) >= this->_minDepth)) {
    795:           this->_size = this->_sieve->support(p)->size();
    796:         } else {
    797:           this->_size = 0;
    798:         }
    799:         return &this->_size;
    800:       };
    801:     public:
    802:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    803:         ostringstream txt;
    804:         int rank;

    806:         if (comm == MPI_COMM_NULL) {
    807:           comm = this->comm();
    808:           rank = this->commRank();
    809:         } else {
    810:           MPI_Comm_rank(comm, &rank);
    811:         }
    812:         if (name == "") {
    813:           if(rank == 0) {
    814:             txt << "viewing a SupportSizeSection" << std::endl;
    815:           }
    816:         } else {
    817:           if(rank == 0) {
    818:             txt << "viewing SupportSizeSection '" << name << "'" << std::endl;
    819:           }
    820:         }
    821:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    822:         PetscSynchronizedFlush(comm);
    823:       };
    824:     };

    826:     template<typename Sieve_>
    827:     class SupportSection : public ALE::ParallelObject {
    828:     public:
    829:       typedef Sieve_                           sieve_type;
    830:       typedef typename sieve_type::source_type point_type;
    831:       typedef typename sieve_type::target_type value_type;
    832:       typedef PartitionDomain<sieve_type>      chart_type;
    833:     protected:
    834:       Obj<sieve_type> _sieve;
    835:       int             _supportSize;
    836:       value_type     *_support;
    837:       chart_type      _domain;
    838:       void ensureSupport(const int size) {
    839:         if (size > this->_supportSize) {
    840:           if (this->_support) delete [] this->_support;
    841:           this->_supportSize = size;
    842:           this->_support     = new value_type[this->_supportSize];
    843:         }
    844:       };
    845:     public:
    846:       SupportSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _supportSize(-1), _support(NULL) {};
    847:       virtual ~SupportSection() {if (this->_support) delete [] this->_support;};
    848:     public:
    849:       const chart_type& getChart() {return this->_domain;};
    850:       bool hasPoint(const point_type& point) {return true;};
    851:       const value_type *restrictPoint(const point_type& p) {
    852:         const Obj<typename sieve_type::traits::supportSequence>& support = this->_sieve->support(p);
    853:         int s = 0;

    855:         this->ensureSupport(support->size());
    856:         for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
    857:           this->_support[s++] = *s_iter;
    858:         }
    859:         return this->_support;
    860:       };
    861:     public:
    862:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    863:         ostringstream txt;
    864:         int rank;

    866:         if (comm == MPI_COMM_NULL) {
    867:           comm = this->comm();
    868:           rank = this->commRank();
    869:         } else {
    870:           MPI_Comm_rank(comm, &rank);
    871:         }
    872:         if (name == "") {
    873:           if(rank == 0) {
    874:             txt << "viewing a SupportSection" << std::endl;
    875:           }
    876:         } else {
    877:           if(rank == 0) {
    878:             txt << "viewing SupportSection '" << name << "'" << std::endl;
    879:           }
    880:         }
    881:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    882:         PetscSynchronizedFlush(comm);
    883:       };
    884:     };

    886:     template<typename Sieve_, typename Section_>
    887:     class ArrowSection : public ALE::ParallelObject {
    888:     public:
    889:       typedef Sieve_                            sieve_type;
    890:       typedef Section_                          section_type;
    891:       typedef typename sieve_type::target_type  point_type;
    892:       typedef typename section_type::point_type arrow_type;
    893:       typedef typename section_type::value_type value_type;
    894:     protected:
    895:       Obj<sieve_type>   _sieve;
    896:       Obj<section_type> _section;
    897:       int               _coneSize;
    898:       value_type       *_cone;
    899:       void ensureCone(const int size) {
    900:         if (size > this->_coneSize) {
    901:           if (this->_cone) delete [] this->_cone;
    902:           this->_coneSize = size;
    903:           this->_cone     = new value_type[this->_coneSize];
    904:         }
    905:       };
    906:     public:
    907:       ArrowSection(const Obj<sieve_type>& sieve, const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _section(section), _coneSize(-1), _cone(NULL) {};
    908:       virtual ~ArrowSection() {if (this->_cone) delete [] this->_cone;};
    909:     public:
    910:       bool hasPoint(const point_type& point) {return this->_sieve->baseContains(point);};
    911:       const value_type *restrictPoint(const point_type& p) {
    912:         const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
    913:         int c = -1;

    915:         this->ensureCone(cone->size());
    916:         for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
    917:           this->_cone[++c] = this->_section->restrictPoint(arrow_type(*c_iter, p))[0];
    918:         }
    919:         return this->_cone;
    920:       };
    921:     public:
    922:       void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
    923:         ostringstream txt;
    924:         int rank;

    926:         if (comm == MPI_COMM_NULL) {
    927:           comm = this->comm();
    928:           rank = this->commRank();
    929:         } else {
    930:           MPI_Comm_rank(comm, &rank);
    931:         }
    932:         if (name == "") {
    933:           if(rank == 0) {
    934:             txt << "viewing a ConeSection" << std::endl;
    935:           }
    936:         } else {
    937:           if(rank == 0) {
    938:             txt << "viewing ConeSection '" << name << "'" << std::endl;
    939:           }
    940:         }
    941:         PetscSynchronizedPrintf(comm, txt.str().c_str());
    942:         PetscSynchronizedFlush(comm);
    943:       };
    944:     };
    945:   }
    946: }
    947: #endif