Actual source code: Sections.hh

  1: #ifndef included_ALE_Sections_hh
  2: #define included_ALE_Sections_hh

  4: #ifndef  included_ALE_Numbering_hh
  5: #include <Numbering.hh>
  6: #endif

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

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

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

110:       if (size > maxSize) {
111:         const value_type dummy(0);

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

140:       for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
141:         array[c++] = *c_iter;
142:       }
143:       return array;
144:     };
145:   };

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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