Actual source code: ISieve.hh
1: #ifndef included_ALE_ISieve_hh
2: #define included_ALE_ISieve_hh
4: #ifndef included_ALE_hh
5: #include <ALE.hh>
6: #endif
8: //#define IMESH_NEW_LABELS
10: namespace ALE {
11: template<typename Point>
12: class OrientedPoint : public std::pair<Point, int> {
13: public:
14: OrientedPoint(const int o) : std::pair<Point, int>(o, o) {};
15: ~OrientedPoint() {};
16: friend std::ostream& operator<<(std::ostream& stream, const OrientedPoint& point) {
17: stream << "(" << point.first << ", " << point.second << ")";
18: return stream;
19: };
20: };
22: template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
23: class Interval {
24: public:
25: typedef Point_ point_type;
26: typedef Alloc_ alloc_type;
27: public:
28: class const_iterator {
29: protected:
30: point_type _p;
31: public:
32: const_iterator(const point_type p): _p(p) {};
33: ~const_iterator() {};
34: public:
35: const_iterator& operator=(const const_iterator& iter) {this->_p = iter._p; return *this;};
36: bool operator==(const const_iterator& iter) const {return this->_p == iter._p;};
37: bool operator!=(const const_iterator& iter) const {return this->_p != iter._p;};
38: const_iterator& operator++() {++this->_p; return *this;}
39: const_iterator& operator++(int) {
40: const_iterator tmp(*this);
41: ++(*this);
42: return tmp;
43: };
44: const_iterator& operator--() {--this->_p; return *this;}
45: const_iterator& operator--(int) {
46: const_iterator tmp(*this);
47: --(*this);
48: return tmp;
49: };
50: point_type operator*() const {return this->_p;};
51: };
52: protected:
53: point_type _min, _max;
54: public:
55: Interval(): _min(point_type()), _max(point_type()) {};
56: Interval(const point_type& min, const point_type& max): _min(min), _max(max) {};
57: Interval(const Interval& interval): _min(interval.min()), _max(interval.max()) {};
58: template<typename Iterator>
59: Interval(Iterator& iterator) {
60: this->_min = *std::min_element(iterator.begin(), iterator.end());
61: this->_max = (*std::max_element(iterator.begin(), iterator.end()))+1;
62: };
63: public:
64: Interval& operator=(const Interval& interval) {_min = interval.min(); _max = interval.max(); return *this;};
65: friend std::ostream& operator<<(std::ostream& stream, const Interval& interval) {
66: stream << "(" << interval.min() << ", " << interval.max() << ")";
67: return stream;
68: };
69: public:
70: const_iterator begin() const {return const_iterator(this->_min);};
71: const_iterator end() const {return const_iterator(this->_max);};
72: size_t size() const {return this->_max - this->_min;};
73: size_t count(const point_type& p) const {return ((p >= _min) && (p < _max)) ? 1 : 0;};
74: point_type min() const {return this->_min;};
75: point_type max() const {return this->_max;};
76: bool hasPoint(const point_type& point) const {
77: if (point < this->_min || point >= this->_max) return false;
78: return true;
79: };
80: void checkPoint(const point_type& point) const {
81: if (point < this->_min || point >= this->_max) {
82: ostringstream msg;
83: msg << "Invalid point " << point << " not in " << *this << std::endl;
84: throw ALE::Exception(msg.str().c_str());
85: }
86: };
87: };
89: template<typename Source_, typename Target_>
90: struct SimpleArrow {
91: typedef Source_ source_type;
92: typedef Target_ target_type;
93: const source_type source;
94: const target_type target;
95: SimpleArrow(const source_type& s, const target_type& t) : source(s), target(t) {};
96: template<typename OtherSource_, typename OtherTarget_>
97: struct rebind {
98: typedef SimpleArrow<OtherSource_, OtherTarget_> other;
99: };
100: struct flip {
101: typedef SimpleArrow<target_type, source_type> other;
102: other arrow(const SimpleArrow& a) {return type(a.target, a.source);};
103: };
104: friend bool operator<(const SimpleArrow& x, const SimpleArrow& y) {
105: return ((x.source < y.source) || ((x.source == y.source) && (x.target < y.target)));
106: };
107: friend std::ostream& operator<<(std::ostream& os, const SimpleArrow& a) {
108: os << a.source << " ----> " << a.target;
109: return os;
110: }
111: };
113: namespace ISieveVisitor {
114: template<typename Sieve>
115: class NullVisitor {
116: public:
117: void visitArrow(const typename Sieve::arrow_type&) {};
118: void visitPoint(const typename Sieve::point_type&) {};
119: void visitArrow(const typename Sieve::arrow_type&, const int orientation) {};
120: void visitPoint(const typename Sieve::point_type&, const int orientation) {};
121: };
122: class PrintVisitor {
123: protected:
124: ostringstream& os;
125: const int rank;
126: public:
127: PrintVisitor(ostringstream& s, const int rank = 0) : os(s), rank(rank) {};
128: template<typename Arrow>
129: void visitArrow(const Arrow& arrow) const {
130: this->os << "["<<this->rank<<"]: " << arrow << std::endl;
131: };
132: template<typename Point>
133: void visitPoint(const Point&) const {};
134: };
135: class ReversePrintVisitor : public PrintVisitor {
136: public:
137: ReversePrintVisitor(ostringstream& s, const int rank) : PrintVisitor(s, rank) {};
138: template<typename Arrow>
139: void visitArrow(const Arrow& arrow) const {
140: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << std::endl;
141: };
142: template<typename Arrow>
143: void visitArrow(const Arrow& arrow, const int orientation) const {
144: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << ": " << orientation << std::endl;
145: };
146: template<typename Point>
147: void visitPoint(const Point&) const {};
148: template<typename Point>
149: void visitPoint(const Point&, const int) const {};
150: };
151: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
152: class PointRetriever {
153: public:
154: typedef typename Sieve::point_type point_type;
155: typedef typename Sieve::arrow_type arrow_type;
156: typedef std::pair<point_type,int> oriented_point_type;
157: protected:
158: const bool unique;
159: size_t i, o;
160: size_t skip, limit;
161: Visitor *visitor;
162: size_t size;
163: point_type *points;
164: oriented_point_type *oPoints;
165: protected:
166: inline virtual bool accept(const point_type& point) {return true;};
167: public:
168: PointRetriever() : unique(false), i(0), o(0), skip(0), limit(0) {
169: this->size = -1;
170: this->points = NULL;
171: this->oPoints = NULL;
172: };
173: PointRetriever(const size_t size, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0) {
174: static Visitor nV;
175: this->visitor = &nV;
176: this->points = NULL;
177: this->oPoints = NULL;
178: this->setSize(size);
179: };
180: PointRetriever(const size_t size, Visitor& v, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0), visitor(&v) {
181: this->points = NULL;
182: this->oPoints = NULL;
183: this->setSize(size);
184: };
185: virtual ~PointRetriever() {
186: delete [] this->points;
187: delete [] this->oPoints;
188: this->points = NULL;
189: this->oPoints = NULL;
190: };
191: void visitArrow(const arrow_type& arrow) {
192: this->visitor->visitArrow(arrow);
193: };
194: void visitArrow(const arrow_type& arrow, const int orientation) {
195: this->visitor->visitArrow(arrow, orientation);
196: };
197: void visitPoint(const point_type& point) {
198: if (i >= size) {
199: ostringstream msg;
200: msg << "Too many points (>" << size << ")for PointRetriever visitor";
201: throw ALE::Exception(msg.str().c_str());
202: }
203: if (this->accept(point)) {
204: if (this->unique) {
205: size_t p;
206: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
207: if (p != i) return;
208: }
209: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
210: points[i++] = point;
211: this->visitor->visitPoint(point);
212: }
213: };
214: void visitPoint(const point_type& point, const int orientation) {
215: if (o >= size) {
216: ostringstream msg;
217: msg << "Too many ordered points (>" << size << ")for PointRetriever visitor";
218: throw ALE::Exception(msg.str().c_str());
219: }
220: if (this->accept(point)) {
221: if (this->unique) {
222: size_t p;
223: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
224: if (p != i) return;
225: }
226: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
227: points[i++] = point;
228: oPoints[o++] = oriented_point_type(point, orientation);
229: this->visitor->visitPoint(point, orientation);
230: }
231: };
232: public:
233: const size_t getSize() const {return this->i;};
234: const point_type *getPoints() const {return this->points;};
235: const size_t getOrientedSize() const {return this->o;};
236: const oriented_point_type *getOrientedPoints() const {return this->oPoints;};
237: void clear() {this->i = this->o = 0;};
238: void setSize(const size_t s) {
239: if (this->points) {
240: delete [] this->points;
241: delete [] this->oPoints;
242: }
243: this->size = s;
244: this->points = new point_type[this->size];
245: this->oPoints = new oriented_point_type[this->size];
246: };
247: void setSkip(size_t s) {this->skip = s;};
248: void setLimit(size_t l) {this->limit = l;};
249: };
250: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
251: class NConeRetriever : public PointRetriever<Sieve,Visitor> {
252: public:
253: typedef PointRetriever<Sieve,Visitor> base_type;
254: typedef typename Sieve::point_type point_type;
255: typedef typename Sieve::arrow_type arrow_type;
256: typedef typename base_type::oriented_point_type oriented_point_type;
257: protected:
258: const Sieve& sieve;
259: protected:
260: inline virtual bool accept(const point_type& point) {
261: if (!this->sieve.getConeSize(point))
262: return true;
263: return false;
264: };
265: public:
266: NConeRetriever(const Sieve& s, const size_t size) : PointRetriever<Sieve,Visitor>(size, true), sieve(s) {};
267: NConeRetriever(const Sieve& s, const size_t size, Visitor& v) : PointRetriever<Sieve,Visitor>(size, v, true), sieve(s) {};
268: virtual ~NConeRetriever() {};
269: };
270: template<typename Mesh, typename Visitor = NullVisitor<typename Mesh::sieve_type> >
271: class MeshNConeRetriever : public PointRetriever<typename Mesh::sieve_type,Visitor> {
272: public:
273: typedef typename Mesh::Sieve Sieve;
274: typedef PointRetriever<Sieve,Visitor> base_type;
275: typedef typename Sieve::point_type point_type;
276: typedef typename Sieve::arrow_type arrow_type;
277: typedef typename base_type::oriented_point_type oriented_point_type;
278: protected:
279: const Mesh& mesh;
280: const int depth;
281: protected:
282: inline virtual bool accept(const point_type& point) {
283: if (this->mesh.depth(point) == this->depth)
284: return true;
285: return false;
286: };
287: public:
288: MeshNConeRetriever(const Mesh& m, const int depth, const size_t size) : PointRetriever<typename Mesh::Sieve,Visitor>(size, true), mesh(m), depth(depth) {};
289: MeshNConeRetriever(const Mesh& m, const int depth, const size_t size, Visitor& v) : PointRetriever<typename Mesh::Sieve,Visitor>(size, v, true), mesh(m), depth(depth) {};
290: virtual ~MeshNConeRetriever() {};
291: };
292: template<typename Sieve, typename Set, typename Renumbering>
293: class FilteredPointRetriever {
294: public:
295: typedef typename Sieve::point_type point_type;
296: typedef typename Sieve::arrow_type arrow_type;
297: typedef std::pair<point_type,int> oriented_point_type;
298: protected:
299: const Set& pointSet;
300: Renumbering& renumbering;
301: const size_t size;
302: size_t i, o;
303: bool renumber;
304: point_type *points;
305: oriented_point_type *oPoints;
306: public:
307: FilteredPointRetriever(const Set& s, Renumbering& r, const size_t size) : pointSet(s), renumbering(r), size(size), i(0), o(0), renumber(true) {
308: this->points = new point_type[this->size];
309: this->oPoints = new oriented_point_type[this->size];
310: };
311: ~FilteredPointRetriever() {
312: delete [] this->points;
313: delete [] this->oPoints;
314: };
315: void visitArrow(const arrow_type& arrow) {};
316: void visitPoint(const point_type& point) {
317: if (i >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
318: if (this->pointSet.find(point) == this->pointSet.end()) return;
319: if (renumber) {
320: points[i++] = this->renumbering[point];
321: } else {
322: points[i++] = point;
323: }
324: };
325: void visitArrow(const arrow_type& arrow, const int orientation) {};
326: void visitPoint(const point_type& point, const int orientation) {
327: if (o >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
328: if (this->pointSet.find(point) == this->pointSet.end()) return;
329: if (renumber) {
330: points[i++] = this->renumbering[point];
331: oPoints[o++] = oriented_point_type(this->renumbering[point], orientation);
332: } else {
333: points[i++] = point;
334: oPoints[o++] = oriented_point_type(point, orientation);
335: }
336: };
337: public:
338: const size_t getSize() const {return this->i;};
339: const point_type *getPoints() const {return this->points;};
340: const size_t getOrientedSize() const {return this->o;};
341: const oriented_point_type *getOrientedPoints() const {return this->oPoints;};
342: void clear() {this->i = 0; this->o = 0;};
343: void useRenumbering(const bool renumber) {this->renumber = renumber;};
344: };
345: template<typename Sieve, int size, typename Visitor = NullVisitor<Sieve> >
346: class ArrowRetriever {
347: public:
348: typedef typename Sieve::point_type point_type;
349: typedef typename Sieve::arrow_type arrow_type;
350: typedef std::pair<arrow_type,int> oriented_arrow_type;
351: protected:
352: arrow_type arrows[size];
353: oriented_arrow_type oArrows[size];
354: size_t i, o;
355: Visitor *visitor;
356: public:
357: ArrowRetriever() : i(0), o(0) {
358: static Visitor nV;
359: this->visitor = &nV;
360: };
361: ArrowRetriever(Visitor& v) : i(0), o(0), visitor(&v) {};
362: void visitArrow(const typename Sieve::arrow_type& arrow) {
363: if (i >= size) throw ALE::Exception("Too many arrows for visitor");
364: arrows[i++] = arrow;
365: this->visitor->visitArrow(arrow);
366: };
367: void visitArrow(const typename Sieve::arrow_type& arrow, const int orientation) {
368: if (o >= size) throw ALE::Exception("Too many arrows for visitor");
369: oArrows[o++] = oriented_arrow_type(arrow, orientation);
370: this->visitor->visitArrow(arrow, orientation);
371: };
372: void visitPoint(const point_type& point) {
373: this->visitor->visitPoint(point);
374: };
375: void visitPoint(const point_type& point, const int orientation) {
376: this->visitor->visitPoint(point, orientation);
377: };
378: public:
379: const size_t getSize() const {return this->i;};
380: const point_type *getArrows() const {return this->arrows;};
381: const size_t getOrientedSize() const {return this->o;};
382: const oriented_arrow_type *getOrientedArrows() const {return this->oArrows;};
383: void clear() {this->i = this->o = 0;};
384: };
385: template<typename Sieve, typename Visitor>
386: class ConeVisitor {
387: protected:
388: const Sieve& sieve;
389: Visitor& visitor;
390: bool useSource;
391: public:
392: ConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
393: void visitPoint(const typename Sieve::point_type& point) {
394: this->sieve.cone(point, visitor);
395: };
396: void visitArrow(const typename Sieve::arrow_type& arrow) {};
397: };
398: template<typename Sieve, typename Visitor>
399: class OrientedConeVisitor {
400: protected:
401: const Sieve& sieve;
402: Visitor& visitor;
403: bool useSource;
404: public:
405: OrientedConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
406: void visitPoint(const typename Sieve::point_type& point) {
407: this->sieve.orientedCone(point, visitor);
408: };
409: void visitArrow(const typename Sieve::arrow_type& arrow) {};
410: };
411: template<typename Sieve, typename Visitor>
412: class SupportVisitor {
413: protected:
414: const Sieve& sieve;
415: Visitor& visitor;
416: bool useSource;
417: public:
418: SupportVisitor(const Sieve& s, Visitor& v, bool useSource = true) : sieve(s), visitor(v), useSource(useSource) {};
419: void visitPoint(const typename Sieve::point_type& point) {
420: this->sieve.support(point, visitor);
421: };
422: void visitArrow(const typename Sieve::arrow_type& arrow) {};
423: };
424: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
425: class TransitiveClosureVisitor {
426: public:
427: typedef Visitor visitor_type;
428: protected:
429: const Sieve& sieve;
430: Visitor& visitor;
431: bool isCone;
432: std::set<typename Sieve::point_type> seen;
433: public:
434: TransitiveClosureVisitor(const Sieve& s, Visitor& v) : sieve(s), visitor(v), isCone(true) {};
435: void visitPoint(const typename Sieve::point_type& point) const {};
436: void visitArrow(const typename Sieve::arrow_type& arrow) {
437: if (this->isCone) {
438: if (this->seen.find(arrow.target) == this->seen.end()) {
439: this->seen.insert(arrow.target);
440: this->visitor.visitPoint(arrow.target);
441: }
442: this->visitor.visitArrow(arrow);
443: if (this->seen.find(arrow.source) == this->seen.end()) {
444: if (this->sieve.getConeSize(arrow.source)) {
445: this->sieve.cone(arrow.source, *this);
446: } else {
447: this->seen.insert(arrow.source);
448: this->visitor.visitPoint(arrow.source);
449: }
450: }
451: } else {
452: if (this->seen.find(arrow.source) == this->seen.end()) {
453: this->seen.insert(arrow.source);
454: this->visitor.visitPoint(arrow.source);
455: }
456: this->visitor.visitArrow(arrow);
457: if (this->seen.find(arrow.target) == this->seen.end()) {
458: if (this->sieve.getSupportSize(arrow.target)) {
459: this->sieve.support(arrow.target, *this);
460: } else {
461: this->seen.insert(arrow.target);
462: this->visitor.visitPoint(arrow.target);
463: }
464: }
465: }
466: };
467: public:
468: bool getIsCone() const {return this->isCone;};
469: void setIsCone(const bool isCone) {this->isCone = isCone;};
470: const std::set<typename Sieve::point_type>& getPoints() const {return this->seen;};
471: void clear() {this->seen.clear();};
472: };
473: template<typename Sieve, typename Section>
474: class SizeVisitor {
475: protected:
476: const Section& section;
477: int size;
478: public:
479: SizeVisitor(const Section& s) : section(s), size(0) {};
480: void visitPoint(const typename Sieve::point_type& point) {
481: this->size += section.getConstrainedFiberDimension(point);
482: };
483: void visitArrow(const typename Sieve::arrow_type&) {};
484: public:
485: int getSize() {return this->size;};
486: };
487: template<typename Sieve, typename Section>
488: class SizeWithBCVisitor {
489: protected:
490: const Section& section;
491: int size;
492: public:
493: SizeWithBCVisitor(const Section& s) : section(s), size(0) {};
494: void visitPoint(const typename Sieve::point_type& point) {
495: this->size += section.getFiberDimension(point);
496: };
497: void visitArrow(const typename Sieve::arrow_type&) {};
498: public:
499: int getSize() {return this->size;};
500: };
501: template<typename Section>
502: class RestrictVisitor {
503: public:
504: typedef typename Section::value_type value_type;
505: protected:
506: const Section& section;
507: int size;
508: int i;
509: value_type *values;
510: bool allocated;
511: public:
512: RestrictVisitor(const Section& s, const int size) : section(s), size(size), i(0) {
513: this->values = new value_type[this->size];
514: this->allocated = true;
515: };
516: RestrictVisitor(const Section& s, const int size, value_type *values) : section(s), size(size), i(0) {
517: this->values = values;
518: this->allocated = false;
519: };
520: ~RestrictVisitor() {if (this->allocated) {delete [] this->values;}};
521: template<typename Point>
522: void visitPoint(const Point& point, const int orientation) {
523: const int dim = section.getFiberDimension(point);
524: if (i+dim > size) {throw ALE::Exception("Too many values for RestrictVisitor.");}
525: const value_type *v = section.restrictPoint(point);
527: if (orientation >= 0) {
528: for(int d = 0; d < dim; ++d, ++i) {
529: this->values[i] = v[d];
530: }
531: } else {
532: for(int d = dim-1; d >= 0; --d, ++i) {
533: this->values[i] = v[d];
534: }
535: }
536: };
537: template<typename Arrow>
538: void visitArrow(const Arrow& arrow, const int orientation) {};
539: public:
540: const value_type *getValues() const {return this->values;};
541: int getSize() const {return this->i;};
542: int getMaxSize() const {return this->size;};
543: void ensureSize(const int size) {
544: this->clear();
545: if (size > this->size) {
546: this->size = size;
547: if (this->allocated) {delete [] this->values;}
548: this->values = new value_type[this->size];
549: this->allocated = true;
550: }
551: };
552: void clear() {this->i = 0;};
553: };
554: template<typename Section>
555: class UpdateVisitor {
556: public:
557: typedef typename Section::value_type value_type;
558: protected:
559: Section& section;
560: const value_type *values;
561: int i;
562: public:
563: UpdateVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
564: template<typename Point>
565: void visitPoint(const Point& point, const int orientation) {
566: const int dim = section.getFiberDimension(point);
567: this->section.updatePoint(point, &this->values[this->i], orientation);
568: this->i += dim;
569: };
570: template<typename Arrow>
571: void visitArrow(const Arrow& arrow, const int orientation) {};
572: };
573: template<typename Section>
574: class UpdateAllVisitor {
575: public:
576: typedef typename Section::value_type value_type;
577: protected:
578: Section& section;
579: const value_type *values;
580: int i;
581: public:
582: UpdateAllVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
583: template<typename Point>
584: void visitPoint(const Point& point, const int orientation) {
585: const int dim = section.getFiberDimension(point);
586: this->section.updatePointAll(point, &this->values[this->i], orientation);
587: this->i += dim;
588: };
589: template<typename Arrow>
590: void visitArrow(const Arrow& arrow, const int orientation) {};
591: };
592: template<typename Section>
593: class UpdateAddVisitor {
594: public:
595: typedef typename Section::value_type value_type;
596: protected:
597: Section& section;
598: const value_type *values;
599: int i;
600: public:
601: UpdateAddVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
602: template<typename Point>
603: void visitPoint(const Point& point, const int orientation) {
604: const int dim = section.getFiberDimension(point);
605: this->section.updateAddPoint(point, &this->values[this->i], orientation);
606: this->i += dim;
607: };
608: template<typename Arrow>
609: void visitArrow(const Arrow& arrow, const int orientation) {};
610: };
611: template<typename Section, typename Order, typename Value>
612: class IndicesVisitor {
613: public:
614: typedef Value value_type;
615: typedef typename Section::point_type point_type;
616: protected:
617: const Section& section;
618: // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
619: Order& order;
620: int size;
621: int i, p;
622: // If false, constrained indices are returned as negative values. Otherwise, they are omitted
623: bool freeOnly;
624: // If true, it allows space for constrained variables (even if the indices are not returned) Wierd
625: bool skipConstraints;
626: value_type *values;
627: bool allocated;
628: point_type *points;
629: bool allocatedPoints;
630: protected:
631: void getUnconstrainedIndices(const point_type& p, const int orientation) {
632: if (i+section.getFiberDimension(p) > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
633: if (orientation >= 0) {
634: const int start = this->order.getIndex(p);
635: const int end = start + section.getFiberDimension(p);
637: for(int j = start; j < end; ++j, ++i) {
638: this->values[i] = j;
639: }
640: } else if (!section.getNumSpaces()) {
641: const int start = this->order.getIndex(p);
642: const int end = start + section.getFiberDimension(p);
644: for(int j = end-1; j >= start; --j, ++i) {
645: this->values[i] = j;
646: }
647: } else {
648: const int numSpaces = section.getNumSpaces();
649: int start = this->order.getIndex(p);
651: for(int space = 0; space < numSpaces; ++space) {
652: const int end = start + section.getFiberDimension(p, space);
654: for(int j = end-1; j >= start; --j, ++i) {
655: this->values[i] = j;
656: }
657: start = end;
658: }
659: }
660: };
661: void getConstrainedIndices(const point_type& p, const int orientation) {
662: const int cDim = this->section.getConstraintDimension(p);
663: if (i+cDim > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
664: typedef typename Section::bc_type::value_type index_type;
665: const index_type *cDof = this->section.getConstraintDof(p);
666: const int start = this->order.getIndex(p);
668: if (orientation >= 0) {
669: const int dim = this->section.getFiberDimension(p);
671: for(int j = start, cInd = 0, k = 0; k < dim; ++k) {
672: if ((cInd < cDim) && (k == cDof[cInd])) {
673: if (!freeOnly) values[i++] = -(k+1);
674: if (skipConstraints) ++j;
675: ++cInd;
676: } else {
677: values[i++] = j++;
678: }
679: }
680: } else {
681: int offset = 0;
682: int cOffset = 0;
683: int k = -1;
685: for(int space = 0; space < section.getNumSpaces(); ++space) {
686: const int dim = this->section.getFiberDimension(p, space);
687: const int tDim = this->section.getConstrainedFiberDimension(p, space);
688: int cInd = (dim - tDim)-1;
690: k += dim;
691: for(int l = 0, j = start+tDim+offset; l < dim; ++l, --k) {
692: if ((cInd >= 0) && (k == cDof[cInd+cOffset])) {
693: if (!freeOnly) values[i++] = -(offset+l+1);
694: if (skipConstraints) --j;
695: --cInd;
696: } else {
697: values[i++] = --j;
698: }
699: }
700: k += dim;
701: offset += dim;
702: cOffset += dim - tDim;
703: }
704: }
705: };
706: public:
707: IndicesVisitor(const Section& s, Order& o, const int size, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
708: this->values = new value_type[this->size];
709: this->allocated = true;
710: if (unique) {
711: this->points = new point_type[this->size];
712: this->allocatedPoints = true;
713: } else {
714: this->points = NULL;
715: this->allocatedPoints = false;
716: }
717: };
718: IndicesVisitor(const Section& s, Order& o, const int size, value_type *values, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
719: this->values = values;
720: this->allocated = false;
721: if (unique) {
722: this->points = new point_type[this->size];
723: this->allocatedPoints = true;
724: } else {
725: this->points = NULL;
726: this->allocatedPoints = false;
727: }
728: };
729: ~IndicesVisitor() {
730: if (this->allocated) {delete [] this->values;}
731: if (this->allocatedPoints) {delete [] this->points;}
732: };
733: void visitPoint(const point_type& point, const int orientation) {
734: if (p >= size) {
735: ostringstream msg;
736: msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
737: throw ALE::Exception(msg.str().c_str());
738: }
739: if (points) {
740: int pp;
741: for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
742: if (pp != p) return;
743: points[p++] = point;
744: }
745: const int cDim = this->section.getConstraintDimension(point);
747: if (!cDim) {
748: this->getUnconstrainedIndices(point, orientation);
749: } else {
750: this->getConstrainedIndices(point, orientation);
751: }
752: };
753: template<typename Arrow>
754: void visitArrow(const Arrow& arrow, const int orientation) {};
755: public:
756: const value_type *getValues() const {return this->values;};
757: int getSize() const {return this->i;};
758: int getMaxSize() const {return this->size;};
759: void ensureSize(const int size) {
760: this->clear();
761: if (size > this->size) {
762: this->size = size;
763: if (this->allocated) {delete [] this->values;}
764: this->values = new value_type[this->size];
765: this->allocated = true;
766: if (this->allocatedPoints) {delete [] this->points;}
767: this->points = new point_type[this->size];
768: this->allocatedPoints = true;
769: }
770: };
771: void clear() {this->i = 0; this->p = 0;};
772: };
773: template<typename Sieve, typename Label>
774: class MarkVisitor {
775: protected:
776: Label& label;
777: int marker;
778: public:
779: MarkVisitor(Label& l, const int marker) : label(l), marker(marker) {};
780: void visitPoint(const typename Sieve::point_type& point) {
781: this->label.setCone(this->marker, point);
782: };
783: void visitArrow(const typename Sieve::arrow_type&) {};
784: };
785: };
787: template<typename Sieve>
788: class ISieveTraversal {
789: public:
790: typedef typename Sieve::point_type point_type;
791: public:
792: template<typename Visitor>
793: static void orientedClosure(const Sieve& sieve, const point_type& p, Visitor& v) {
794: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
795: Retriever cV[2] = {Retriever(200,v), Retriever(200,v)};
796: int c = 0;
798: v.visitPoint(p, 0);
799: // Cone is guarateed to be ordered correctly
800: sieve.orientedCone(p, cV[c]);
802: while(cV[c].getOrientedSize()) {
803: const typename Retriever::oriented_point_type *cone = cV[c].getOrientedPoints();
804: const int coneSize = cV[c].getOrientedSize();
805: c = 1 - c;
807: for(int p = 0; p < coneSize; ++p) {
808: const typename Retriever::point_type& point = cone[p].first;
809: int pO = cone[p].second == 0 ? 1 : cone[p].second;
810: const int pConeSize = sieve.getConeSize(point);
812: if (pO < 0) {
813: if (pO == -pConeSize) {
814: sieve.orientedReverseCone(point, cV[c]);
815: } else {
816: const int numSkip = sieve.getConeSize(point) + pO;
818: cV[c].setSkip(cV[c].getSize()+numSkip);
819: cV[c].setLimit(cV[c].getSize()+pConeSize);
820: sieve.orientedReverseCone(point, cV[c]);
821: sieve.orientedReverseCone(point, cV[c]);
822: cV[c].setSkip(0);
823: cV[c].setLimit(0);
824: }
825: } else {
826: if (pO == 1) {
827: sieve.orientedCone(point, cV[c]);
828: } else {
829: const int numSkip = pO-1;
831: cV[c].setSkip(cV[c].getSize()+numSkip);
832: cV[c].setLimit(cV[c].getSize()+pConeSize);
833: sieve.orientedCone(point, cV[c]);
834: sieve.orientedCone(point, cV[c]);
835: cV[c].setSkip(0);
836: cV[c].setLimit(0);
837: }
838: }
839: }
840: cV[1-c].clear();
841: }
842: #if 0
843: // These contain arrows paired with orientations from the \emph{previous} arrow
844: Obj<orientedArrowArray> cone = new orientedArrowArray();
845: Obj<orientedArrowArray> base = new orientedArrowArray();
846: coneSet seen;
848: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
849: const arrow_type arrow(*c_iter, p);
851: cone->push_back(oriented_arrow_type(arrow, 1)); // Notice the orientation of leaf cells is always 1
852: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0])); // However, we return the actual arrow orientation
853: }
854: for(int i = 1; i < depth; ++i) {
855: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
857: cone->clear();
858: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
859: const arrow_type& arrow = b_iter->first; // This is an arrow considered in the previous round
860: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source); // We are going to get the cone of this guy
861: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0]; // The orientation of arrow, which is our pO
863: if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
864: o = -(o+1);
865: }
866: if (o < 0) {
867: const int size = pCone->size();
869: if (o == -size) {
870: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
871: if (seen.find(*c_iter) == seen.end()) {
872: const arrow_type newArrow(*c_iter, arrow.source);
873: int pointO = orientation->restrictPoint(newArrow)[0];
875: seen.insert(*c_iter);
876: cone->push_back(oriented_arrow_type(newArrow, o));
877: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
878: }
879: }
880: } else {
881: const int numSkip = size + o;
882: int count = 0;
884: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
885: if (count < numSkip) continue;
886: if (seen.find(*c_iter) == seen.end()) {
887: const arrow_type newArrow(*c_iter, arrow.source);
888: int pointO = orientation->restrictPoint(newArrow)[0];
890: seen.insert(*c_iter);
891: cone->push_back(oriented_arrow_type(newArrow, o));
892: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
893: }
894: }
895: count = 0;
896: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
897: if (count >= numSkip) break;
898: if (seen.find(*c_iter) == seen.end()) {
899: const arrow_type newArrow(*c_iter, arrow.source);
900: int pointO = orientation->restrictPoint(newArrow)[0];
902: seen.insert(*c_iter);
903: cone->push_back(oriented_arrow_type(newArrow, o));
904: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
905: }
906: }
907: }
908: } else {
909: if (o == 1) {
910: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
911: if (seen.find(*c_iter) == seen.end()) {
912: const arrow_type newArrow(*c_iter, arrow.source);
914: seen.insert(*c_iter);
915: cone->push_back(oriented_arrow_type(newArrow, o));
916: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
917: }
918: }
919: } else {
920: const int numSkip = o-1;
921: int count = 0;
923: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
924: if (count < numSkip) continue;
925: if (seen.find(*c_iter) == seen.end()) {
926: const arrow_type newArrow(*c_iter, arrow.source);
928: seen.insert(*c_iter);
929: cone->push_back(oriented_arrow_type(newArrow, o));
930: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
931: }
932: }
933: count = 0;
934: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
935: if (count >= numSkip) break;
936: if (seen.find(*c_iter) == seen.end()) {
937: const arrow_type newArrow(*c_iter, arrow.source);
939: seen.insert(*c_iter);
940: cone->push_back(oriented_arrow_type(newArrow, o));
941: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
942: }
943: }
944: }
945: }
946: }
947: }
948: #endif
949: };
950: };
952: namespace IFSieveDef {
953: template<typename PointType_>
954: class Sequence {
955: public:
956: typedef PointType_ point_type;
957: class const_iterator {
958: public:
959: // Standard iterator typedefs
960: typedef std::input_iterator_tag iterator_category;
961: typedef PointType_ value_type;
962: typedef int difference_type;
963: typedef value_type* pointer;
964: typedef value_type& reference;
965: protected:
966: const point_type *_data;
967: int _pos;
968: public:
969: const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
970: ~const_iterator() {};
971: public:
972: virtual bool operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
973: virtual bool operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
974: virtual const value_type operator*() const {return this->_data[this->_pos];};
975: virtual const_iterator& operator++() {++this->_pos; return *this;};
976: virtual const_iterator operator++(int n) {
977: const_iterator tmp(*this);
978: ++this->_pos;
979: return tmp;
980: };
981: };
982: typedef const_iterator iterator;
983: protected:
984: const point_type *_data;
985: int _begin;
986: int _end;
987: public:
988: Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
989: ~Sequence() {};
990: public:
991: virtual iterator begin() const {return iterator(_data, _begin);};
992: virtual iterator end() const {return iterator(_data, _end);};
993: size_t size() const {return(_end - _begin);}
994: void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
995: };
996: }
998: // Interval Final Sieve
999: // This is just two CSR matrices that give cones and supports
1000: // It is completely static and cannot be resized
1001: // It will operator on visitors, rather than sequences (which are messy)
1002: template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1003: class IFSieve : public ParallelObject {
1004: public:
1005: // Types
1006: typedef IFSieve<Point_,Allocator_> this_type;
1007: typedef Point_ point_type;
1008: typedef SimpleArrow<point_type,point_type> arrow_type;
1009: typedef typename arrow_type::source_type source_type;
1010: typedef typename arrow_type::target_type target_type;
1011: typedef int index_type;
1012: // Allocators
1013: typedef Allocator_ point_allocator_type;
1014: typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1015: typedef typename point_allocator_type::template rebind<int>::other int_allocator_type;
1016: // Interval
1017: typedef Interval<point_type, point_allocator_type> chart_type;
1018: // Dynamic structure
1019: typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1020: // Iterator interface
1021: typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1022: typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1023: // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1024: typedef std::set<point_type> pointSet;
1025: typedef ALE::array<point_type> pointArray;
1026: typedef pointSet coneSet;
1027: typedef pointSet supportSet;
1028: typedef pointArray coneArray;
1029: typedef pointArray supportArray;
1030: protected:
1031: // Arrow Containers
1032: typedef index_type *offsets_type;
1033: typedef point_type *cones_type;
1034: typedef point_type *supports_type;
1035: // Decorators
1036: typedef int *orientations_type;
1037: protected:
1038: // Data
1039: bool indexAllocated;
1040: offsets_type coneOffsets;
1041: offsets_type supportOffsets;
1042: bool pointAllocated;
1043: index_type maxConeSize;
1044: index_type maxSupportSize;
1045: index_type baseSize;
1046: index_type capSize;
1047: cones_type cones;
1048: supports_type supports;
1049: bool orientCones;
1050: orientations_type coneOrientations;
1051: chart_type chart;
1052: int_allocator_type intAlloc;
1053: index_allocator_type indexAlloc;
1054: point_allocator_type pointAlloc;
1055: newpoints_type newCones;
1056: newpoints_type newSupports;
1057: // Sequences
1058: Obj<coneSequence> coneSeq;
1059: Obj<supportSequence> supportSeq;
1060: protected: // Memory Management
1061: void createIndices() {
1062: this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1063: this->coneOffsets -= this->chart.min();
1064: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1065: this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1066: this->supportOffsets -= this->chart.min();
1067: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1068: this->indexAllocated = true;
1069: };
1070: void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1071: if (*coneOffsets) {
1072: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1073: *coneOffsets += chart.min();
1074: indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1075: *coneOffsets = NULL;
1076: }
1077: if (*supportOffsets) {
1078: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1079: *supportOffsets += chart.min();
1080: indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1081: *supportOffsets = NULL;
1082: }
1083: };
1084: void destroyIndices() {
1085: this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1086: this->indexAllocated = false;
1087: this->maxConeSize = -1;
1088: this->maxSupportSize = -1;
1089: this->baseSize = -1;
1090: this->capSize = -1;
1091: };
1092: void createPoints() {
1093: this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1094: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1095: this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1096: for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1097: if (orientCones) {
1098: this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1099: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1100: }
1101: this->pointAllocated = true;
1102: };
1103: void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1104: if (*cones) {
1105: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1106: pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1107: *cones = NULL;
1108: }
1109: if (*supports) {
1110: for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1111: pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1112: *supports = NULL;
1113: }
1114: if (*coneOrientations) {
1115: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1116: intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1117: *coneOrientations = NULL;
1118: }
1119: };
1120: void destroyPoints() {
1121: this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1122: this->pointAllocated = false;
1123: };
1124: void prefixSum(const offsets_type array) {
1125: for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1126: array[p] = array[p] + array[p-1];
1127: }
1128: };
1129: void calculateBaseAndCapSize() {
1130: this->baseSize = 0;
1131: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1132: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1133: ++this->baseSize;
1134: }
1135: }
1136: this->capSize = 0;
1137: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1138: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1139: ++this->capSize;
1140: }
1141: }
1142: };
1143: public:
1144: IFSieve(const MPI_Comm comm, const int debug = 0) : ParallelObject(comm, debug), indexAllocated(false), coneOffsets(NULL), supportOffsets(NULL), pointAllocated(false), maxConeSize(-1), maxSupportSize(-1), baseSize(-1), capSize(-1), cones(NULL), supports(NULL), orientCones(true), coneOrientations(NULL) {
1145: this->coneSeq = new coneSequence(NULL, 0, 0);
1146: this->supportSeq = new supportSequence(NULL, 0, 0);
1147: };
1148: IFSieve(const MPI_Comm comm, const point_type& min, const point_type& max, const int debug = 0) : ParallelObject(comm, debug), indexAllocated(false), coneOffsets(NULL), supportOffsets(NULL), pointAllocated(false), maxConeSize(-1), maxSupportSize(-1), baseSize(-1), capSize(-1), cones(NULL), supports(NULL), orientCones(true), coneOrientations(NULL) {
1149: this->setChart(chart_type(min, max));
1150: this->coneSeq = new coneSequence(NULL, 0, 0);
1151: this->supportSeq = new supportSequence(NULL, 0, 0);
1152: };
1153: ~IFSieve() {
1154: this->destroyPoints();
1155: this->destroyIndices();
1156: };
1157: public: // Accessors
1158: const chart_type& getChart() const {return this->chart;};
1159: void setChart(const chart_type& chart) {
1160: this->destroyPoints();
1161: this->destroyIndices();
1162: this->chart = chart;
1163: this->createIndices();
1164: };
1165: index_type getMaxConeSize() const {
1166: return this->maxConeSize;
1167: };
1168: index_type getMaxSupportSize() const {
1169: return this->maxSupportSize;
1170: };
1171: bool baseContains(const point_type& p) const {
1172: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1173: this->chart.checkPoint(p);
1175: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1176: return true;
1177: }
1178: return false;
1179: };
1180: public: // Construction
1181: index_type getConeSize(const point_type& p) const {
1182: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1183: this->chart.checkPoint(p);
1184: return this->coneOffsets[p+1]-this->coneOffsets[p];
1185: };
1186: void setConeSize(const point_type& p, const index_type c) {
1187: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1188: this->chart.checkPoint(p);
1189: this->coneOffsets[p+1] = c;
1190: this->maxConeSize = std::max(this->maxConeSize, c);
1191: };
1192: void addConeSize(const point_type& p, const index_type c) {
1193: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1194: this->chart.checkPoint(p);
1195: this->coneOffsets[p+1] += c;
1196: this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1197: };
1198: index_type getSupportSize(const point_type& p) const {
1199: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1200: this->chart.checkPoint(p);
1201: return this->supportOffsets[p+1]-this->supportOffsets[p];
1202: };
1203: void setSupportSize(const point_type& p, const index_type s) {
1204: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1205: this->chart.checkPoint(p);
1206: this->supportOffsets[p+1] = s;
1207: this->maxSupportSize = std::max(this->maxSupportSize, s);
1208: };
1209: void addSupportSize(const point_type& p, const index_type s) {
1210: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1211: this->chart.checkPoint(p);
1212: this->supportOffsets[p+1] += s;
1213: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1214: };
1215: void allocate() {
1216: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1217: this->prefixSum(this->coneOffsets);
1218: this->prefixSum(this->supportOffsets);
1219: this->createPoints();
1220: this->calculateBaseAndCapSize();
1221: };
1222: // Purely for backwards compatibility
1223: template<typename Color>
1224: void addArrow(const point_type& p, const point_type& q, const Color c) {
1225: this->addArrow(p, q);
1226: };
1227: void addArrow(const point_type& p, const point_type& q) {
1228: if (!this->chart.hasPoint(q)) {
1229: if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1230: const index_type start = this->coneOffsets[q];
1231: const index_type end = this->coneOffsets[q+1];
1233: for(int c = start; c < end; ++c) {
1234: this->newCones[q].push_back(this->cones[c]);
1235: }
1236: }
1237: this->newCones[q].push_back(p);
1238: }
1239: if (!this->chart.hasPoint(p)) {
1240: if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1241: const index_type start = this->supportOffsets[p];
1242: const index_type end = this->supportOffsets[p+1];
1244: for(int s = start; s < end; ++s) {
1245: this->newSupports[p].push_back(this->supports[s]);
1246: }
1247: }
1248: this->newSupports[p].push_back(q);
1249: }
1250: };
1251: void reallocate() {
1252: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1253: if (!this->newCones.size() && !this->newSupports.size()) return;
1254: const chart_type oldChart = this->chart;
1255: offsets_type oldConeOffsets = this->coneOffsets;
1256: offsets_type oldSupportOffsets = this->supportOffsets;
1257: cones_type oldCones = this->cones;
1258: supports_type oldSupports = this->supports;
1259: orientations_type oldConeOrientations = this->coneOrientations;
1260: point_type min = this->chart.min();
1261: point_type max = this->chart.max()-1;
1263: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1264: min = std::min(min, c_iter->first);
1265: max = std::max(max, c_iter->first);
1266: }
1267: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1268: min = std::min(min, s_iter->first);
1269: max = std::max(max, s_iter->first);
1270: }
1271: this->chart = chart_type(min, max+1);
1272: this->createIndices();
1273: // Copy sizes (converted from offsets)
1274: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1275: this->coneOffsets[p+1] = oldConeOffsets[p+1]-oldConeOffsets[p];
1276: this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1277: }
1278: // Inject new sizes
1279: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1280: this->coneOffsets[c_iter->first+1] = c_iter->second.size();
1281: this->maxConeSize = std::max(this->maxConeSize, (int) c_iter->second.size());
1282: }
1283: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1284: this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1285: this->maxSupportSize = std::max(this->maxSupportSize, (int) s_iter->second.size());
1286: }
1287: this->prefixSum(this->coneOffsets);
1288: this->prefixSum(this->supportOffsets);
1289: this->createPoints();
1290: this->calculateBaseAndCapSize();
1291: // Copy cones and supports
1292: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1293: const index_type cStart = this->coneOffsets[p];
1294: const index_type cOStart = oldConeOffsets[p];
1295: const index_type cOEnd = oldConeOffsets[p+1];
1296: const index_type sStart = this->supportOffsets[p];
1297: const index_type sOStart = oldSupportOffsets[p];
1298: const index_type sOEnd = oldSupportOffsets[p+1];
1300: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1301: this->cones[c] = oldCones[cO];
1302: }
1303: for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1304: this->supports[s] = oldSupports[sO];
1305: }
1306: if (this->orientCones) {
1307: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1308: this->coneOrientations[c] = oldConeOrientations[cO];
1309: }
1310: }
1311: }
1312: // Inject new cones and supports
1313: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1314: index_type start = this->coneOffsets[c_iter->first];
1316: for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1317: this->cones[start++] = *p_iter;
1318: }
1319: if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1320: }
1321: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1322: index_type start = this->supportOffsets[s_iter->first];
1324: for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1325: this->supports[start++] = *p_iter;
1326: }
1327: if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1328: }
1329: this->newCones.clear();
1330: this->newSupports.clear();
1331: this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1332: this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1333: };
1334: // Recalculate the support offsets and fill the supports
1335: // This is used when an IFSieve is being used as a label
1336: void recalculateLabel() {
1337: ISieveVisitor::PointRetriever<IFSieve> v(1);
1339: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1340: this->supportOffsets[p+1] = 0;
1341: }
1342: this->maxSupportSize = 0;
1343: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1344: this->cone(p, v);
1345: if (v.getSize()) {
1346: this->supportOffsets[v.getPoints()[0]+1]++;
1347: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1348: }
1349: v.clear();
1350: }
1351: this->prefixSum(this->supportOffsets);
1352: this->calculateBaseAndCapSize();
1353: this->symmetrize();
1354: };
1355: void setCone(const point_type cone[], const point_type& p) {
1356: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1357: this->chart.checkPoint(p);
1358: const index_type start = this->coneOffsets[p];
1359: const index_type end = this->coneOffsets[p+1];
1361: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1362: this->cones[c] = cone[i];
1363: }
1364: };
1365: void setCone(const point_type cone, const point_type& p) {
1366: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1367: this->chart.checkPoint(p);
1368: const index_type start = this->coneOffsets[p];
1369: const index_type end = this->coneOffsets[p+1];
1371: if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1372: this->cones[start] = cone;
1373: };
1374: #if 0
1375: template<typename PointSequence>
1376: void setCone(const PointSequence& cone, const point_type& p) {
1377: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1378: this->chart.checkPoint(p);
1379: const index_type start = this->coneOffsets[p];
1380: const index_type end = this->coneOffsets[p+1];
1381: if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1382: typename PointSequence::iterator c_iter = cone.begin();
1384: for(index_type c = start; c < end; ++c, ++c_iter) {
1385: this->cones[c] = *c_iter;
1386: }
1387: };
1388: #endif
1389: void setSupport(const point_type& p, const point_type support[]) {
1390: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1391: this->chart.checkPoint(p);
1392: const index_type start = this->supportOffsets[p];
1393: const index_type end = this->supportOffsets[p+1];
1395: for(index_type s = start, i = 0; s < end; ++s, ++i) {
1396: this->supports[s] = support[i];
1397: }
1398: };
1399: #if 0
1400: template<typename PointSequence>
1401: void setSupport(const point_type& p, const PointSequence& support) {
1402: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1403: this->chart.checkPoint(p);
1404: const index_type start = this->supportOffsets[p];
1405: const index_type end = this->supportOffsets[p+1];
1406: if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1407: typename PointSequence::iterator s_iter = support.begin();
1409: for(index_type s = start; s < end; ++s, ++s_iter) {
1410: this->supports[s] = *s_iter;
1411: }
1412: };
1413: #endif
1414: void setConeOrientation(const int coneOrientation[], const point_type& p) {
1415: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1416: this->chart.checkPoint(p);
1417: const index_type start = this->coneOffsets[p];
1418: const index_type end = this->coneOffsets[p+1];
1420: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1421: this->coneOrientations[c] = coneOrientation[i];
1422: }
1423: };
1424: void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1425: for(point_type p = 0; p < numCells; ++p) {
1426: const index_type start = p*numCorners;
1427: const index_type end = (p+1)*numCorners;
1429: for(index_type c = start; c < end; ++c) {
1430: const point_type q = cones[c]+offset;
1432: this->supportOffsets[q+1]++;
1433: }
1434: }
1435: for(point_type p = numCells; p < this->chart.max(); ++p) {
1436: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1437: }
1438: };
1439: void symmetrize() {
1440: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1441: offsets -= this->chart.min();
1442: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1443: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1445: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1446: const index_type start = this->coneOffsets[p];
1447: const index_type end = this->coneOffsets[p+1];
1449: for(index_type c = start; c < end; ++c) {
1450: const point_type q = this->cones[c];
1452: this->chart.checkPoint(q);
1453: this->supports[this->supportOffsets[q]+offsets[q]] = p;
1454: ++offsets[q];
1455: }
1456: }
1457: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1458: indexAlloc.deallocate(offsets, this->chart.size()+1);
1459: };
1460: index_type getBaseSize() const {
1461: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1462: return this->baseSize;
1463: };
1464: index_type getCapSize() const {
1465: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1466: return this->capSize;
1467: };
1468: public: // Traversals
1469: template<typename Visitor>
1470: void roots(const Visitor& v) const {
1471: this->roots(const_cast<Visitor&>(v));
1472: };
1473: template<typename Visitor>
1474: void roots(Visitor& v) const {
1475: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1477: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1478: if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1479: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1480: v.visitPoint(p);
1481: }
1482: }
1483: }
1484: };
1485: template<typename Visitor>
1486: void leaves(const Visitor& v) const {
1487: this->leaves(const_cast<Visitor&>(v));
1488: };
1489: template<typename Visitor>
1490: void leaves(Visitor& v) const {
1491: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1493: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1494: if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1495: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1496: v.visitPoint(p);
1497: }
1498: }
1499: }
1500: };
1501: template<typename Visitor>
1502: void base(const Visitor& v) const {
1503: this->base(const_cast<Visitor&>(v));
1504: };
1505: template<typename Visitor>
1506: void base(Visitor& v) const {
1507: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1509: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1510: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1511: v.visitPoint(p);
1512: }
1513: }
1514: };
1515: template<typename Visitor>
1516: void cap(const Visitor& v) const {
1517: this->cap(const_cast<Visitor&>(v));
1518: };
1519: template<typename Visitor>
1520: void cap(Visitor& v) const {
1521: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1523: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1524: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1525: v.visitPoint(p);
1526: }
1527: }
1528: };
1529: template<typename PointSequence, typename Visitor>
1530: void cone(const PointSequence& points, Visitor& v) const {
1531: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1532: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1533: const point_type p = *p_iter;
1534: this->chart.checkPoint(p);
1535: const index_type start = this->coneOffsets[p];
1536: const index_type end = this->coneOffsets[p+1];
1538: for(index_type c = start; c < end; ++c) {
1539: v.visitArrow(arrow_type(this->cones[c], p));
1540: v.visitPoint(this->cones[c]);
1541: }
1542: }
1543: };
1544: template<typename Visitor>
1545: void cone(const point_type& p, Visitor& v) const {
1546: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1547: this->chart.checkPoint(p);
1548: const index_type start = this->coneOffsets[p];
1549: const index_type end = this->coneOffsets[p+1];
1551: for(index_type c = start; c < end; ++c) {
1552: v.visitArrow(arrow_type(this->cones[c], p));
1553: v.visitPoint(this->cones[c]);
1554: }
1555: };
1556: const Obj<coneSequence>& cone(const point_type& p) const {
1557: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1558: if (!this->chart.hasPoint(p)) {
1559: this->coneSeq->reset(this->cones, 0, 0);
1560: } else {
1561: this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1562: }
1563: return this->coneSeq;
1564: };
1565: template<typename PointSequence, typename Visitor>
1566: void support(const PointSequence& points, Visitor& v) const {
1567: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1568: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1569: const point_type p = *p_iter;
1570: this->chart.checkPoint(p);
1571: const index_type start = this->supportOffsets[p];
1572: const index_type end = this->supportOffsets[p+1];
1574: for(index_type s = start; s < end; ++s) {
1575: v.visitArrow(arrow_type(p, this->supports[s]));
1576: v.visitPoint(this->supports[s]);
1577: }
1578: }
1579: };
1580: template<typename Visitor>
1581: void support(const point_type& p, Visitor& v) const {
1582: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1583: this->chart.checkPoint(p);
1584: const index_type start = this->supportOffsets[p];
1585: const index_type end = this->supportOffsets[p+1];
1587: for(index_type s = start; s < end; ++s) {
1588: v.visitArrow(arrow_type(p, this->supports[s]));
1589: v.visitPoint(this->supports[s]);
1590: }
1591: };
1592: const Obj<supportSequence>& support(const point_type& p) const {
1593: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1594: if (!this->chart.hasPoint(p)) {
1595: this->supportSeq->reset(this->supports, 0, 0);
1596: } else {
1597: this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1598: }
1599: return this->supportSeq;
1600: };
1601: template<typename Visitor>
1602: void orientedCone(const point_type& p, Visitor& v) const {
1603: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1604: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1605: this->chart.checkPoint(p);
1606: const index_type start = this->coneOffsets[p];
1607: const index_type end = this->coneOffsets[p+1];
1609: for(index_type c = start; c < end; ++c) {
1610: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1611: v.visitPoint(this->cones[c], this->coneOrientations[c]);
1612: }
1613: };
1614: template<typename Visitor>
1615: void orientedReverseCone(const point_type& p, Visitor& v) const {
1616: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1617: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1618: this->chart.checkPoint(p);
1619: const index_type start = this->coneOffsets[p];
1620: const index_type end = this->coneOffsets[p+1];
1622: for(index_type c = end-1; c >= start; --c) {
1623: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1624: v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1625: }
1626: };
1627: // Currently does only 1 level
1628: // Does not check for uniqueness
1629: template<typename Visitor>
1630: void meet(const point_type& p, const point_type& q, Visitor& v) const {
1631: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1632: this->chart.checkPoint(p);
1633: this->chart.checkPoint(q);
1634: const index_type startP = this->coneOffsets[p];
1635: const index_type endP = this->coneOffsets[p+1];
1636: const index_type startQ = this->coneOffsets[q];
1637: const index_type endQ = this->coneOffsets[q+1];
1639: for(index_type cP = startP; cP < endP; ++cP) {
1640: const point_type& c1 = this->cones[cP];
1642: for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1643: if (c1 == this->cones[cQ]) v.visitPoint(c1);
1644: }
1645: if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1646: }
1647: };
1648: // Currently does only 1 level
1649: template<typename Sequence, typename Visitor>
1650: void join(const Sequence& points, Visitor& v) const {
1651: typedef std::set<point_type> pointSet;
1652: pointSet intersect[2] = {pointSet(), pointSet()};
1653: pointSet tmp;
1654: int p = 0;
1655: int c = 0;
1657: for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1658: this->chart.checkPoint(*p_iter);
1659: tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1660: if (p == 0) {
1661: intersect[1-c].insert(tmp.begin(), tmp.end());
1662: p++;
1663: } else {
1664: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1665: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1666: intersect[c].clear();
1667: }
1668: c = 1 - c;
1669: tmp.clear();
1670: }
1671: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1672: v.visitPoint(*p_iter);
1673: }
1674: };
1675: public: // Viewing
1676: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1677: ostringstream txt;
1678: int rank;
1680: if (comm == MPI_COMM_NULL) {
1681: comm = this->comm();
1682: rank = this->commRank();
1683: } else {
1684: MPI_Comm_rank(comm, &rank);
1685: }
1686: if (name == "") {
1687: if(rank == 0) {
1688: txt << "viewing an IFSieve" << std::endl;
1689: }
1690: } else {
1691: if(rank == 0) {
1692: txt << "viewing IFSieve '" << name << "'" << std::endl;
1693: }
1694: }
1695: if(rank == 0) {
1696: txt << "cap --> base:" << std::endl;
1697: }
1698: ISieveVisitor::PrintVisitor pV(txt, rank);
1699: this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1700: PetscSynchronizedPrintf(comm, txt.str().c_str());
1701: PetscSynchronizedFlush(comm);
1702: ostringstream txt2;
1704: if(rank == 0) {
1705: txt2 << "base <-- cap:" << std::endl;
1706: }
1707: ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1708: this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1709: PetscSynchronizedPrintf(comm, txt2.str().c_str());
1710: PetscSynchronizedFlush(comm);
1711: if (orientCones) {
1712: ostringstream txt3;
1714: if(rank == 0) {
1715: txt3 << "Orientation:" << std::endl;
1716: }
1717: ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1718: this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1719: PetscSynchronizedPrintf(comm, txt3.str().c_str());
1720: PetscSynchronizedFlush(comm);
1721: }
1722: };
1723: };
1725: class ISieveConverter {
1726: public:
1727: template<typename Sieve, typename ISieve, typename Renumbering>
1728: static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1729: // First construct a renumbering of the sieve points
1730: const Obj<typename Sieve::baseSequence>& base = sieve.base();
1731: const Obj<typename Sieve::capSequence>& cap = sieve.cap();
1732: typename ISieve::point_type min = 0;
1733: typename ISieve::point_type max = 0;
1735: if (renumber) {
1736: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1737: renumbering[*b_iter] = max++;
1738: }
1739: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1740: if (renumbering.find(*c_iter) == renumbering.end()) {
1741: renumbering[*c_iter] = max++;
1742: }
1743: }
1744: } else {
1745: if (base->size()) {
1746: min = *base->begin();
1747: max = *base->begin();
1748: } else if (cap->size()) {
1749: min = *cap->begin();
1750: max = *cap->begin();
1751: }
1752: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1753: min = std::min(min, *b_iter);
1754: max = std::max(max, *b_iter);
1755: }
1756: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1757: min = std::min(min, *c_iter);
1758: max = std::max(max, *c_iter);
1759: }
1760: if (base->size() || cap->size()) {
1761: ++max;
1762: }
1763: for(typename ISieve::point_type p = min; p < max; ++p) {
1764: renumbering[p] = p;
1765: }
1766: }
1767: // Create the ISieve
1768: isieve.setChart(typename ISieve::chart_type(min, max));
1769: // Set cone and support sizes
1770: size_t maxSize = 0;
1772: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1773: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1775: isieve.setConeSize(renumbering[*b_iter], cone->size());
1776: maxSize = std::max(maxSize, cone->size());
1777: }
1778: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1779: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1781: isieve.setSupportSize(renumbering[*c_iter], support->size());
1782: maxSize = std::max(maxSize, support->size());
1783: }
1784: isieve.allocate();
1785: // Fill up cones and supports
1786: typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];
1788: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1789: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1790: int i = 0;
1792: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1793: points[i] = renumbering[*c_iter];
1794: }
1795: isieve.setCone(points, renumbering[*b_iter]);
1796: }
1797: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1798: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1799: int i = 0;
1801: for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
1802: points[i] = renumbering[*s_iter];
1803: }
1804: isieve.setSupport(renumbering[*c_iter], points);
1805: }
1806: delete [] points;
1807: };
1808: template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
1809: static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
1810: if (isieve.getMaxConeSize() < 0) return;
1811: const Obj<typename Sieve::baseSequence>& base = sieve.base();
1812: int *orientations = new int[isieve.getMaxConeSize()];
1814: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1815: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1816: int i = 0;
1818: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1819: typename ArrowSection::point_type arrow(*c_iter, *b_iter);
1821: orientations[i] = orientation->restrictPoint(arrow)[0];
1822: }
1823: isieve.setConeOrientation(orientations, renumbering[*b_iter]);
1824: }
1825: delete [] orientations;
1826: };
1827: template<typename Section, typename ISection, typename Renumbering>
1828: static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
1829: const typename Section::chart_type& chart = coordinates.getChart();
1830: typename ISection::point_type min = *chart.begin();
1831: typename ISection::point_type max = *chart.begin();
1833: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1834: min = std::min(min, *p_iter);
1835: max = std::max(max, *p_iter);
1836: }
1837: icoordinates.setChart(typename ISection::chart_type(min, max+1));
1838: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1839: icoordinates.setFiberDimension(*p_iter, coordinates.getFiberDimension(*p_iter));
1840: }
1841: icoordinates.allocatePoint();
1842: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1843: icoordinates.updatePoint(*p_iter, coordinates.restrictPoint(*p_iter));
1844: }
1845: };
1846: template<typename IMesh, typename Label>
1847: static void convertLabel(IMesh& imesh, const std::string& name, const Obj<Label>& oldLabel) {
1848: const Obj<typename IMesh::label_type>& label = imesh.createLabel(name);
1849: const typename IMesh::sieve_type::chart_type& chart = imesh.getSieve()->getChart();
1850: int size = 0;
1852: label->setChart(chart);
1853: for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1854: const int coneSize = oldLabel->cone(p)->size();
1856: label->setConeSize(p, coneSize);
1857: size += coneSize;
1858: }
1859: if (size) {label->setSupportSize(0, size);}
1860: label->allocate();
1861: for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1862: const Obj<typename Label::coneSequence>& cone = oldLabel->cone(p);
1864: if (cone->size()) {
1865: label->setCone(*cone->begin(), p);
1866: }
1867: }
1868: label->recalculateLabel();
1869: };
1870: template<typename Mesh, typename IMesh, typename Renumbering>
1871: static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
1872: convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
1873: imesh.stratify();
1874: convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
1875: convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
1876: const typename Mesh::labels_type& labels = mesh.getLabels();
1877: std::string heightName("height");
1878: std::string depthName("depth");
1880: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
1881: #ifdef IMESH_NEW_LABELS
1882: if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
1883: convertLabel(imesh, l_iter->first, l_iter->second);
1884: }
1885: #else
1886: imesh.setLabel(l_iter->first, l_iter->second);
1887: #endif
1888: }
1889: };
1890: };
1891: }
1893: #endif