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: #include <fstream>
10: //#define IMESH_NEW_LABELS
12: namespace ALE {
13: template<typename Point>
14: class OrientedPoint : public std::pair<Point, int> {
15: public:
16: OrientedPoint(const int o) : std::pair<Point, int>(o, o) {};
17: ~OrientedPoint() {};
18: friend std::ostream& operator<<(std::ostream& stream, const OrientedPoint& point) {
19: stream << "(" << point.first << ", " << point.second << ")";
20: return stream;
21: };
22: };
24: template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
25: class Interval {
26: public:
27: typedef Point_ point_type;
28: typedef Alloc_ alloc_type;
29: public:
30: class const_iterator {
31: protected:
32: point_type _p;
33: public:
34: const_iterator(const point_type p): _p(p) {};
35: ~const_iterator() {};
36: public:
37: const_iterator& operator=(const const_iterator& iter) {this->_p = iter._p; return *this;};
38: bool operator==(const const_iterator& iter) const {return this->_p == iter._p;};
39: bool operator!=(const const_iterator& iter) const {return this->_p != iter._p;};
40: const_iterator& operator++() {++this->_p; return *this;}
41: const_iterator& operator++(int) {
42: const_iterator tmp(*this);
43: ++(*this);
44: return tmp;
45: };
46: const_iterator& operator--() {--this->_p; return *this;}
47: const_iterator& operator--(int) {
48: const_iterator tmp(*this);
49: --(*this);
50: return tmp;
51: };
52: point_type operator*() const {return this->_p;};
53: };
54: protected:
55: point_type _min, _max;
56: public:
57: Interval(): _min(point_type()), _max(point_type()) {};
58: Interval(const point_type& min, const point_type& max): _min(min), _max(max) {};
59: Interval(const Interval& interval): _min(interval.min()), _max(interval.max()) {};
60: template<typename Iterator>
61: Interval(Iterator& iterator) {
62: this->_min = *std::min_element(iterator.begin(), iterator.end());
63: this->_max = (*std::max_element(iterator.begin(), iterator.end()))+1;
64: }
65: public:
66: Interval& operator=(const Interval& interval) {_min = interval.min(); _max = interval.max(); return *this;}
67: friend std::ostream& operator<<(std::ostream& stream, const Interval& interval) {
68: stream << "(" << interval.min() << ", " << interval.max() << ")";
69: return stream;
70: }
71: public:
72: const_iterator begin() const {return const_iterator(this->_min);};
73: const_iterator end() const {return const_iterator(this->_max);};
74: size_t size() const {return this->_max - this->_min;};
75: size_t count(const point_type& p) const {return ((p >= _min) && (p < _max)) ? 1 : 0;};
76: point_type min() const {return this->_min;};
77: point_type max() const {return this->_max;};
78: bool hasPoint(const point_type& point) const {
79: if (point < this->_min || point >= this->_max) return false;
80: return true;
81: };
82: void checkPoint(const point_type& point) const {
83: if (point < this->_min || point >= this->_max) {
84: ostringstream msg;
85: msg << "Invalid point " << point << " not in " << *this << std::endl;
86: throw ALE::Exception(msg.str().c_str());
87: }
88: };
89: };
91: template<typename Source_, typename Target_>
92: struct SimpleArrow {
93: typedef Source_ source_type;
94: typedef Target_ target_type;
95: const source_type source;
96: const target_type target;
97: SimpleArrow(const source_type& s, const target_type& t) : source(s), target(t) {};
98: template<typename OtherSource_, typename OtherTarget_>
99: struct rebind {
100: typedef SimpleArrow<OtherSource_, OtherTarget_> other;
101: };
102: struct flip {
103: typedef SimpleArrow<target_type, source_type> other;
104: other arrow(const SimpleArrow& a) {return type(a.target, a.source);};
105: };
106: friend bool operator<(const SimpleArrow& x, const SimpleArrow& y) {
107: return ((x.source < y.source) || ((x.source == y.source) && (x.target < y.target)));
108: };
109: friend std::ostream& operator<<(std::ostream& os, const SimpleArrow& a) {
110: os << a.source << " ----> " << a.target;
111: return os;
112: }
113: };
115: namespace ISieveVisitor {
116: template<typename Sieve>
117: class NullVisitor {
118: public:
119: inline void visitArrow(const typename Sieve::arrow_type&) {};
120: inline void visitPoint(const typename Sieve::point_type&) {};
121: inline void visitArrow(const typename Sieve::arrow_type&, const int orientation) {};
122: inline void visitPoint(const typename Sieve::point_type&, const int orientation) {};
123: };
124: class PrintVisitor {
125: protected:
126: ostringstream& os;
127: const int rank;
128: public:
129: PrintVisitor(ostringstream& s, const int rank = 0) : os(s), rank(rank) {};
130: template<typename Arrow>
131: inline void visitArrow(const Arrow& arrow) const {
132: this->os << "["<<this->rank<<"]: " << arrow << std::endl;
133: }
134: template<typename Point>
135: inline void visitPoint(const Point&) const {}
136: };
137: class ReversePrintVisitor : public PrintVisitor {
138: public:
139: ReversePrintVisitor(ostringstream& s, const int rank) : PrintVisitor(s, rank) {};
140: template<typename Arrow>
141: inline void visitArrow(const Arrow& arrow) const {
142: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << std::endl;
143: }
144: template<typename Arrow>
145: inline void visitArrow(const Arrow& arrow, const int orientation) const {
146: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << ": " << orientation << std::endl;
147: }
148: template<typename Point>
149: inline void visitPoint(const Point&) const {}
150: template<typename Point>
151: inline void visitPoint(const Point&, const int) const {}
152: };
153: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
154: class PointRetriever {
155: public:
156: typedef typename Sieve::point_type point_type;
157: typedef typename Sieve::arrow_type arrow_type;
158: typedef std::pair<point_type,int> oriented_point_type;
159: protected:
160: const bool unique;
161: size_t i, o;
162: size_t skip, limit;
163: Visitor *visitor;
164: size_t size;
165: point_type *points;
166: oriented_point_type *oPoints;
167: protected:
168: inline virtual bool accept(const point_type& point) {return true;};
169: public:
170: PointRetriever() : unique(false), i(0), o(0), skip(0), limit(0) {
171: this->size = 0;
172: this->points = NULL;
173: this->oPoints = NULL;
174: };
175: PointRetriever(const size_t size, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0) {
176: static Visitor nV;
177: this->visitor = &nV;
178: this->points = NULL;
179: this->oPoints = NULL;
180: this->setSize(size);
181: };
182: PointRetriever(const size_t size, Visitor& v, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0), visitor(&v) {
183: this->points = NULL;
184: this->oPoints = NULL;
185: this->setSize(size);
186: };
187: virtual ~PointRetriever() {
188: delete [] this->points;
189: delete [] this->oPoints;
190: this->points = NULL;
191: this->oPoints = NULL;
192: };
193: inline void visitArrow(const arrow_type& arrow) {
194: this->visitor->visitArrow(arrow);
195: };
196: inline void visitArrow(const arrow_type& arrow, const int orientation) {
197: this->visitor->visitArrow(arrow, orientation);
198: };
199: inline void visitPoint(const point_type& point) {
200: if (i >= size) {
201: ostringstream msg;
202: msg << "Too many points (>" << size << ")for PointRetriever visitor";
203: throw ALE::Exception(msg.str().c_str());
204: }
205: if (this->accept(point)) {
206: if (this->unique) {
207: size_t p;
208: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
209: if (p != i) return;
210: }
211: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
212: points[i++] = point;
213: this->visitor->visitPoint(point);
214: }
215: };
216: inline void visitPoint(const point_type& point, const int orientation) {
217: if (o >= size) {
218: ostringstream msg;
219: msg << "Too many ordered points (>" << size << ")for PointRetriever visitor";
220: throw ALE::Exception(msg.str().c_str());
221: }
222: if (this->accept(point)) {
223: if (this->unique) {
224: size_t p;
225: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
226: if (p != i) return;
227: }
228: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
229: points[i++] = point;
230: oPoints[o++] = oriented_point_type(point, orientation);
231: this->visitor->visitPoint(point, orientation);
232: }
233: };
234: public:
235: size_t getSize() const {return this->i;}
236: const point_type *getPoints() const {return this->points;}
237: size_t getOrientedSize() const {return this->o;}
238: const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
239: void clear() {this->i = this->o = 0;}
240: void setSize(const size_t s) {
241: if (this->points) {
242: delete [] this->points;
243: delete [] this->oPoints;
244: }
245: this->size = s;
246: this->points = new point_type[this->size];
247: this->oPoints = new oriented_point_type[this->size];
248: }
249: void setSkip(size_t s) {this->skip = s;};
250: void setLimit(size_t l) {this->limit = l;};
251: };
252: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
253: class NConeRetriever : public PointRetriever<Sieve,Visitor> {
254: public:
255: typedef PointRetriever<Sieve,Visitor> base_type;
256: typedef typename Sieve::point_type point_type;
257: typedef typename Sieve::arrow_type arrow_type;
258: typedef typename base_type::oriented_point_type oriented_point_type;
259: protected:
260: const Sieve& sieve;
261: protected:
262: inline virtual bool accept(const point_type& point) {
263: if (!this->sieve.getConeSize(point))
264: return true;
265: return false;
266: };
267: public:
268: NConeRetriever(const Sieve& s, const size_t size) : PointRetriever<Sieve,Visitor>(size, true), sieve(s) {};
269: NConeRetriever(const Sieve& s, const size_t size, Visitor& v) : PointRetriever<Sieve,Visitor>(size, v, true), sieve(s) {};
270: virtual ~NConeRetriever() {};
271: };
272: template<typename Mesh, typename Visitor = NullVisitor<typename Mesh::sieve_type> >
273: class MeshNConeRetriever : public PointRetriever<typename Mesh::sieve_type,Visitor> {
274: public:
275: typedef typename Mesh::Sieve Sieve;
276: typedef PointRetriever<Sieve,Visitor> base_type;
277: typedef typename Sieve::point_type point_type;
278: typedef typename Sieve::arrow_type arrow_type;
279: typedef typename base_type::oriented_point_type oriented_point_type;
280: protected:
281: const Mesh& mesh;
282: const int depth;
283: protected:
284: inline virtual bool accept(const point_type& point) {
285: if (this->mesh.depth(point) == this->depth)
286: return true;
287: return false;
288: };
289: public:
290: MeshNConeRetriever(const Mesh& m, const int depth, const size_t size) : PointRetriever<typename Mesh::Sieve,Visitor>(size, true), mesh(m), depth(depth) {};
291: MeshNConeRetriever(const Mesh& m, const int depth, const size_t size, Visitor& v) : PointRetriever<typename Mesh::Sieve,Visitor>(size, v, true), mesh(m), depth(depth) {};
292: virtual ~MeshNConeRetriever() {};
293: };
294: template<typename Sieve, typename Set, typename Renumbering>
295: class FilteredPointRetriever {
296: public:
297: typedef typename Sieve::point_type point_type;
298: typedef typename Sieve::arrow_type arrow_type;
299: typedef std::pair<point_type,int> oriented_point_type;
300: protected:
301: const Set& pointSet;
302: Renumbering& renumbering;
303: const size_t size;
304: size_t i, o;
305: bool renumber;
306: point_type *points;
307: oriented_point_type *oPoints;
308: public:
309: FilteredPointRetriever(const Set& s, Renumbering& r, const size_t size) : pointSet(s), renumbering(r), size(size), i(0), o(0), renumber(true) {
310: this->points = new point_type[this->size];
311: this->oPoints = new oriented_point_type[this->size];
312: };
313: ~FilteredPointRetriever() {
314: delete [] this->points;
315: delete [] this->oPoints;
316: };
317: inline void visitArrow(const arrow_type& arrow) {};
318: inline void visitPoint(const point_type& point) {
319: if (i >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
320: if (this->pointSet.find(point) == this->pointSet.end()) return;
321: if (renumber) {
322: points[i++] = this->renumbering[point];
323: } else {
324: points[i++] = point;
325: }
326: };
327: inline void visitArrow(const arrow_type& arrow, const int orientation) {};
328: inline void visitPoint(const point_type& point, const int orientation) {
329: if (o >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
330: if (this->pointSet.find(point) == this->pointSet.end()) return;
331: if (renumber) {
332: points[i++] = this->renumbering[point];
333: oPoints[o++] = oriented_point_type(this->renumbering[point], orientation);
334: } else {
335: points[i++] = point;
336: oPoints[o++] = oriented_point_type(point, orientation);
337: }
338: };
339: public:
340: size_t getSize() const {return this->i;}
341: const point_type *getPoints() const {return this->points;}
342: size_t getOrientedSize() const {return this->o;}
343: const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
344: void clear() {this->i = 0; this->o = 0;}
345: void useRenumbering(const bool renumber) {this->renumber = renumber;}
346: };
347: template<typename Sieve, int size, typename Visitor = NullVisitor<Sieve> >
348: class ArrowRetriever {
349: public:
350: typedef typename Sieve::point_type point_type;
351: typedef typename Sieve::arrow_type arrow_type;
352: typedef std::pair<arrow_type,int> oriented_arrow_type;
353: protected:
354: arrow_type arrows[size];
355: oriented_arrow_type oArrows[size];
356: size_t i, o;
357: Visitor *visitor;
358: public:
359: ArrowRetriever() : i(0), o(0) {
360: static Visitor nV;
361: this->visitor = &nV;
362: };
363: ArrowRetriever(Visitor& v) : i(0), o(0), visitor(&v) {};
364: inline void visitArrow(const typename Sieve::arrow_type& arrow) {
365: if (i >= size) throw ALE::Exception("Too many arrows for visitor");
366: arrows[i++] = arrow;
367: this->visitor->visitArrow(arrow);
368: };
369: inline void visitArrow(const typename Sieve::arrow_type& arrow, const int orientation) {
370: if (o >= size) throw ALE::Exception("Too many arrows for visitor");
371: oArrows[o++] = oriented_arrow_type(arrow, orientation);
372: this->visitor->visitArrow(arrow, orientation);
373: };
374: inline void visitPoint(const point_type& point) {
375: this->visitor->visitPoint(point);
376: };
377: inline void visitPoint(const point_type& point, const int orientation) {
378: this->visitor->visitPoint(point, orientation);
379: };
380: public:
381: size_t getSize() const {return this->i;}
382: const point_type *getArrows() const {return this->arrows;}
383: size_t getOrientedSize() const {return this->o;}
384: const oriented_arrow_type *getOrientedArrows() const {return this->oArrows;}
385: void clear() {this->i = this->o = 0;}
386: };
387: template<typename Sieve, typename Visitor>
388: class ConeVisitor {
389: protected:
390: const Sieve& sieve;
391: Visitor& visitor;
392: bool useSource;
393: public:
394: ConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
395: inline void visitPoint(const typename Sieve::point_type& point) {
396: this->sieve.cone(point, visitor);
397: };
398: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
399: };
400: template<typename Sieve, typename Visitor>
401: class OrientedConeVisitor {
402: protected:
403: const Sieve& sieve;
404: Visitor& visitor;
405: bool useSource;
406: public:
407: OrientedConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
408: inline void visitPoint(const typename Sieve::point_type& point) {
409: this->sieve.orientedCone(point, visitor);
410: };
411: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
412: };
413: template<typename Sieve, typename Visitor>
414: class SupportVisitor {
415: protected:
416: const Sieve& sieve;
417: Visitor& visitor;
418: bool useSource;
419: public:
420: SupportVisitor(const Sieve& s, Visitor& v, bool useSource = true) : sieve(s), visitor(v), useSource(useSource) {};
421: inline void visitPoint(const typename Sieve::point_type& point) {
422: this->sieve.support(point, visitor);
423: };
424: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
425: };
426: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
427: class TransitiveClosureVisitor {
428: public:
429: typedef Visitor visitor_type;
430: protected:
431: const Sieve& sieve;
432: Visitor& visitor;
433: bool isCone;
434: std::set<typename Sieve::point_type> seen;
435: public:
436: TransitiveClosureVisitor(const Sieve& s, Visitor& v) : sieve(s), visitor(v), isCone(true) {};
437: inline void visitPoint(const typename Sieve::point_type& point) const {};
438: inline void visitArrow(const typename Sieve::arrow_type& arrow) {
439: if (this->isCone) {
440: if (this->seen.find(arrow.target) == this->seen.end()) {
441: this->seen.insert(arrow.target);
442: this->visitor.visitPoint(arrow.target);
443: }
444: this->visitor.visitArrow(arrow);
445: if (this->seen.find(arrow.source) == this->seen.end()) {
446: if (this->sieve.getConeSize(arrow.source)) {
447: this->sieve.cone(arrow.source, *this);
448: } else {
449: this->seen.insert(arrow.source);
450: this->visitor.visitPoint(arrow.source);
451: }
452: }
453: } else {
454: if (this->seen.find(arrow.source) == this->seen.end()) {
455: this->seen.insert(arrow.source);
456: this->visitor.visitPoint(arrow.source);
457: }
458: this->visitor.visitArrow(arrow);
459: if (this->seen.find(arrow.target) == this->seen.end()) {
460: if (this->sieve.getSupportSize(arrow.target)) {
461: this->sieve.support(arrow.target, *this);
462: } else {
463: this->seen.insert(arrow.target);
464: this->visitor.visitPoint(arrow.target);
465: }
466: }
467: }
468: };
469: public:
470: bool getIsCone() const {return this->isCone;};
471: void setIsCone(const bool isCone) {this->isCone = isCone;};
472: const std::set<typename Sieve::point_type>& getPoints() const {return this->seen;};
473: void clear() {this->seen.clear();};
474: };
475: template<typename Sieve, typename Section>
476: class SizeVisitor {
477: protected:
478: const Section& section;
479: int size;
480: public:
481: SizeVisitor(const Section& s) : section(s), size(0) {};
482: inline void visitPoint(const typename Sieve::point_type& point) {
483: this->size += section.getConstrainedFiberDimension(point);
484: };
485: inline void visitArrow(const typename Sieve::arrow_type&) {};
486: public:
487: int getSize() {return this->size;};
488: };
489: template<typename Sieve, typename Section>
490: class SizeWithBCVisitor {
491: protected:
492: const Section& section;
493: int size;
494: public:
495: SizeWithBCVisitor(const Section& s) : section(s), size(0) {};
496: inline void visitPoint(const typename Sieve::point_type& point) {
497: this->size += section.getFiberDimension(point);
498: };
499: inline void visitArrow(const typename Sieve::arrow_type&) {};
500: public:
501: int getSize() {return this->size;};
502: };
503: template<typename Section>
504: class RestrictVisitor {
505: public:
506: typedef typename Section::value_type value_type;
507: protected:
508: const Section& section;
509: int size;
510: int i;
511: value_type *values;
512: bool allocated;
513: public:
514: RestrictVisitor(const Section& s, const int size) : section(s), size(size), i(0) {
515: this->values = new value_type[this->size];
516: this->allocated = true;
517: };
518: RestrictVisitor(const Section& s, const int size, value_type *values) : section(s), size(size), i(0) {
519: this->values = values;
520: this->allocated = false;
521: };
522: ~RestrictVisitor() {if (this->allocated) {delete [] this->values;}};
523: template<typename Point>
524: inline void visitPoint(const Point& point, const int orientation) {
525: const int dim = section.getFiberDimension(point);
526: if (i+dim > size) {throw ALE::Exception("Too many values for RestrictVisitor.");}
527: const value_type *v = section.restrictPoint(point);
529: if (orientation >= 0) {
530: for(int d = 0; d < dim; ++d, ++i) {
531: this->values[i] = v[d];
532: }
533: } else {
534: for(int d = dim-1; d >= 0; --d, ++i) {
535: this->values[i] = v[d];
536: }
537: }
538: }
539: template<typename Arrow>
540: inline void visitArrow(const Arrow& arrow, const int orientation) {}
541: public:
542: const value_type *getValues() const {return this->values;};
543: int getSize() const {return this->i;};
544: int getMaxSize() const {return this->size;};
545: void ensureSize(const int size) {
546: this->clear();
547: if (size > this->size) {
548: this->size = size;
549: if (this->allocated) {delete [] this->values;}
550: this->values = new value_type[this->size];
551: this->allocated = true;
552: }
553: };
554: void clear() {this->i = 0;};
555: };
556: template<typename Section>
557: class UpdateVisitor {
558: public:
559: typedef typename Section::value_type value_type;
560: protected:
561: Section& section;
562: const value_type *values;
563: int i;
564: public:
565: UpdateVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
566: template<typename Point>
567: inline void visitPoint(const Point& point, const int orientation) {
568: const int dim = section.getFiberDimension(point);
569: this->section.updatePoint(point, &this->values[this->i], orientation);
570: this->i += dim;
571: }
572: template<typename Arrow>
573: inline void visitArrow(const Arrow& arrow, const int orientation) {}
574: void clear() {this->i = 0;};
575: };
576: template<typename Section>
577: class UpdateAllVisitor {
578: public:
579: typedef typename Section::value_type value_type;
580: protected:
581: Section& section;
582: const value_type *values;
583: int i;
584: public:
585: UpdateAllVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
586: template<typename Point>
587: inline void visitPoint(const Point& point, const int orientation) {
588: const int dim = section.getFiberDimension(point);
589: this->section.updatePointAll(point, &this->values[this->i], orientation);
590: this->i += dim;
591: }
592: template<typename Arrow>
593: inline void visitArrow(const Arrow& arrow, const int orientation) {}
594: void clear() {this->i = 0;};
595: };
596: template<typename Section>
597: class UpdateAddVisitor {
598: public:
599: typedef typename Section::value_type value_type;
600: protected:
601: Section& section;
602: const value_type *values;
603: int i;
604: public:
605: UpdateAddVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
606: template<typename Point>
607: inline void visitPoint(const Point& point, const int orientation) {
608: const int dim = section.getFiberDimension(point);
609: this->section.updateAddPoint(point, &this->values[this->i], orientation);
610: this->i += dim;
611: }
612: template<typename Arrow>
613: inline void visitArrow(const Arrow& arrow, const int orientation) {}
614: void clear() {this->i = 0;};
615: };
616: template<typename Section, typename Order, typename Value>
617: class IndicesVisitor {
618: public:
619: typedef Value value_type;
620: typedef typename Section::point_type point_type;
621: protected:
622: const Section& section;
623: // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
624: Order& order;
625: int size;
626: int i, p;
627: // If false, constrained indices are returned as negative values. Otherwise, they are omitted
628: bool freeOnly;
629: // If true, it allows space for constrained variables (even if the indices are not returned) Wierd
630: bool skipConstraints;
631: value_type *values;
632: bool allocated;
633: point_type *points;
634: bool allocatedPoints;
635: protected:
636: void getUnconstrainedIndices(const point_type& p, const int orientation) {
637: if (i+section.getFiberDimension(p) > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
638: if (orientation >= 0) {
639: const int start = this->order.getIndex(p);
640: const int end = start + section.getFiberDimension(p);
642: for(int j = start; j < end; ++j, ++i) {
643: this->values[i] = j;
644: }
645: } else if (!section.getNumSpaces()) {
646: const int start = this->order.getIndex(p);
647: const int end = start + section.getFiberDimension(p);
649: for(int j = end-1; j >= start; --j, ++i) {
650: this->values[i] = j;
651: }
652: } else {
653: const int numSpaces = section.getNumSpaces();
654: int start = this->order.getIndex(p);
656: for(int space = 0; space < numSpaces; ++space) {
657: const int end = start + section.getFiberDimension(p, space);
659: for(int j = end-1; j >= start; --j, ++i) {
660: this->values[i] = j;
661: }
662: start = end;
663: }
664: }
665: };
666: void getConstrainedIndices(const point_type& p, const int orientation) {
667: const int cDim = this->section.getConstraintDimension(p);
668: if (i+cDim > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
669: typedef typename Section::bc_type::value_type index_type;
670: const index_type *cDof = this->section.getConstraintDof(p);
671: const int start = this->order.getIndex(p);
673: if (orientation >= 0) {
674: const int dim = this->section.getFiberDimension(p);
676: for(int j = start, cInd = 0, k = 0; k < dim; ++k) {
677: if ((cInd < cDim) && (k == cDof[cInd])) {
678: if (!freeOnly) values[i++] = -(k+1);
679: if (skipConstraints) ++j;
680: ++cInd;
681: } else {
682: values[i++] = j++;
683: }
684: }
685: } else {
686: int offset = 0;
687: int cOffset = 0;
688: int k = -1;
690: for(int space = 0; space < section.getNumSpaces(); ++space) {
691: const int dim = this->section.getFiberDimension(p, space);
692: const int tDim = this->section.getConstrainedFiberDimension(p, space);
693: int cInd = (dim - tDim)-1;
695: k += dim;
696: for(int l = 0, j = start+tDim+offset; l < dim; ++l, --k) {
697: if ((cInd >= 0) && (k == cDof[cInd+cOffset])) {
698: if (!freeOnly) values[i++] = -(offset+l+1);
699: if (skipConstraints) --j;
700: --cInd;
701: } else {
702: values[i++] = --j;
703: }
704: }
705: k += dim;
706: offset += dim;
707: cOffset += dim - tDim;
708: }
709: }
710: };
711: public:
712: IndicesVisitor(const Section& s, Order& o, const int size, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
713: this->values = new value_type[this->size];
714: this->allocated = true;
715: if (unique) {
716: this->points = new point_type[this->size];
717: this->allocatedPoints = true;
718: } else {
719: this->points = NULL;
720: this->allocatedPoints = false;
721: }
722: };
723: IndicesVisitor(const Section& s, Order& o, const int size, value_type *values, const bool unique = false) : section(s), order(o), size(size), i(0), p(0), freeOnly(false), skipConstraints(false) {
724: this->values = values;
725: this->allocated = false;
726: if (unique) {
727: this->points = new point_type[this->size];
728: this->allocatedPoints = true;
729: } else {
730: this->points = NULL;
731: this->allocatedPoints = false;
732: }
733: };
734: ~IndicesVisitor() {
735: if (this->allocated) {delete [] this->values;}
736: if (this->allocatedPoints) {delete [] this->points;}
737: };
738: inline void visitPoint(const point_type& point, const int orientation) {
739: if (p >= size) {
740: ostringstream msg;
741: msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
742: throw ALE::Exception(msg.str().c_str());
743: }
744: if (points) {
745: int pp;
746: for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
747: if (pp != p) return;
748: points[p++] = point;
749: }
750: const int cDim = this->section.getConstraintDimension(point);
752: if (!cDim) {
753: this->getUnconstrainedIndices(point, orientation);
754: } else {
755: this->getConstrainedIndices(point, orientation);
756: }
757: }
758: template<typename Arrow>
759: inline void visitArrow(const Arrow& arrow, const int orientation) {}
760: public:
761: const value_type *getValues() const {return this->values;};
762: int getSize() const {return this->i;};
763: int getMaxSize() const {return this->size;};
764: void ensureSize(const int size) {
765: this->clear();
766: if (size > this->size) {
767: this->size = size;
768: if (this->allocated) {delete [] this->values;}
769: this->values = new value_type[this->size];
770: this->allocated = true;
771: if (this->allocatedPoints) {delete [] this->points;}
772: this->points = new point_type[this->size];
773: this->allocatedPoints = true;
774: }
775: };
776: void clear() {this->i = 0; this->p = 0;};
777: };
778: template<typename Sieve, typename Label>
779: class MarkVisitor {
780: protected:
781: Label& label;
782: int marker;
783: public:
784: MarkVisitor(Label& l, const int marker) : label(l), marker(marker) {};
785: inline void visitPoint(const typename Sieve::point_type& point) {
786: this->label.setCone(this->marker, point);
787: };
788: inline void visitArrow(const typename Sieve::arrow_type&) {};
789: };
790: };
792: template<typename Sieve>
793: class ISieveTraversal {
794: public:
795: typedef typename Sieve::point_type point_type;
796: public:
797: template<typename Visitor>
798: static void orientedClosure(const Sieve& sieve, const point_type& p, Visitor& v) {
799: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
800: Retriever cV[2] = {Retriever(200,v), Retriever(200,v)};
801: int c = 0;
803: v.visitPoint(p, 0);
804: // Cone is guarateed to be ordered correctly
805: sieve.orientedCone(p, cV[c]);
807: while(cV[c].getOrientedSize()) {
808: const typename Retriever::oriented_point_type *cone = cV[c].getOrientedPoints();
809: const int coneSize = cV[c].getOrientedSize();
810: c = 1 - c;
812: for(int p = 0; p < coneSize; ++p) {
813: const typename Retriever::point_type& point = cone[p].first;
814: int pO = cone[p].second == 0 ? 1 : cone[p].second;
815: const int pConeSize = sieve.getConeSize(point);
817: if (pO < 0) {
818: if (pO == -pConeSize) {
819: sieve.orientedReverseCone(point, cV[c]);
820: } else {
821: const int numSkip = sieve.getConeSize(point) + pO;
823: cV[c].setSkip(cV[c].getSize()+numSkip);
824: cV[c].setLimit(cV[c].getSize()+pConeSize);
825: sieve.orientedReverseCone(point, cV[c]);
826: sieve.orientedReverseCone(point, cV[c]);
827: cV[c].setSkip(0);
828: cV[c].setLimit(0);
829: }
830: } else {
831: if (pO == 1) {
832: sieve.orientedCone(point, cV[c]);
833: } else {
834: const int numSkip = pO-1;
836: cV[c].setSkip(cV[c].getSize()+numSkip);
837: cV[c].setLimit(cV[c].getSize()+pConeSize);
838: sieve.orientedCone(point, cV[c]);
839: sieve.orientedCone(point, cV[c]);
840: cV[c].setSkip(0);
841: cV[c].setLimit(0);
842: }
843: }
844: }
845: cV[1-c].clear();
846: }
847: #if 0
848: // These contain arrows paired with orientations from the \emph{previous} arrow
849: Obj<orientedArrowArray> cone = new orientedArrowArray();
850: Obj<orientedArrowArray> base = new orientedArrowArray();
851: coneSet seen;
853: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
854: const arrow_type arrow(*c_iter, p);
856: cone->push_back(oriented_arrow_type(arrow, 1)); // Notice the orientation of leaf cells is always 1
857: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0])); // However, we return the actual arrow orientation
858: }
859: for(int i = 1; i < depth; ++i) {
860: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
862: cone->clear();
863: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
864: const arrow_type& arrow = b_iter->first; // This is an arrow considered in the previous round
865: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source); // We are going to get the cone of this guy
866: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0]; // The orientation of arrow, which is our pO
868: if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
869: o = -(o+1);
870: }
871: if (o < 0) {
872: const int size = pCone->size();
874: if (o == -size) {
875: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
876: if (seen.find(*c_iter) == seen.end()) {
877: const arrow_type newArrow(*c_iter, arrow.source);
878: int pointO = orientation->restrictPoint(newArrow)[0];
880: seen.insert(*c_iter);
881: cone->push_back(oriented_arrow_type(newArrow, o));
882: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
883: }
884: }
885: } else {
886: const int numSkip = size + o;
887: int count = 0;
889: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
890: if (count < numSkip) continue;
891: if (seen.find(*c_iter) == seen.end()) {
892: const arrow_type newArrow(*c_iter, arrow.source);
893: int pointO = orientation->restrictPoint(newArrow)[0];
895: seen.insert(*c_iter);
896: cone->push_back(oriented_arrow_type(newArrow, o));
897: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
898: }
899: }
900: count = 0;
901: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
902: if (count >= numSkip) break;
903: if (seen.find(*c_iter) == seen.end()) {
904: const arrow_type newArrow(*c_iter, arrow.source);
905: int pointO = orientation->restrictPoint(newArrow)[0];
907: seen.insert(*c_iter);
908: cone->push_back(oriented_arrow_type(newArrow, o));
909: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
910: }
911: }
912: }
913: } else {
914: if (o == 1) {
915: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
916: if (seen.find(*c_iter) == seen.end()) {
917: const arrow_type newArrow(*c_iter, arrow.source);
919: seen.insert(*c_iter);
920: cone->push_back(oriented_arrow_type(newArrow, o));
921: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
922: }
923: }
924: } else {
925: const int numSkip = o-1;
926: int count = 0;
928: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
929: if (count < numSkip) continue;
930: if (seen.find(*c_iter) == seen.end()) {
931: const arrow_type newArrow(*c_iter, arrow.source);
933: seen.insert(*c_iter);
934: cone->push_back(oriented_arrow_type(newArrow, o));
935: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
936: }
937: }
938: count = 0;
939: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
940: if (count >= numSkip) break;
941: if (seen.find(*c_iter) == seen.end()) {
942: const arrow_type newArrow(*c_iter, arrow.source);
944: seen.insert(*c_iter);
945: cone->push_back(oriented_arrow_type(newArrow, o));
946: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
947: }
948: }
949: }
950: }
951: }
952: }
953: #endif
954: }
955: template<typename Visitor>
956: static void orientedStar(const Sieve& sieve, const point_type& p, Visitor& v) {
957: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
958: Retriever sV[2] = {Retriever(200,v), Retriever(200,v)};
959: int s = 0;
961: v.visitPoint(p, 0);
962: // Support is guarateed to be ordered correctly
963: sieve.orientedSupport(p, sV[s]);
965: while(sV[s].getOrientedSize()) {
966: const typename Retriever::oriented_point_type *support = sV[s].getOrientedPoints();
967: const int supportSize = sV[s].getOrientedSize();
968: s = 1 - s;
970: for(int p = 0; p < supportSize; ++p) {
971: const typename Retriever::point_type& point = support[p].first;
972: int pO = support[p].second == 0 ? 1 : support[p].second;
973: const int pSupportSize = sieve.getSupportSize(point);
975: if (pO < 0) {
976: if (pO == -pSupportSize) {
977: sieve.orientedReverseSupport(point, sV[s]);
978: } else {
979: const int numSkip = sieve.getSupportSize(point) + pO;
981: sV[s].setSkip(sV[s].getSize()+numSkip);
982: sV[s].setLimit(sV[s].getSize()+pSupportSize);
983: sieve.orientedReverseSupport(point, sV[s]);
984: sieve.orientedReverseSupport(point, sV[s]);
985: sV[s].setSkip(0);
986: sV[s].setLimit(0);
987: }
988: } else {
989: if (pO == 1) {
990: sieve.orientedSupport(point, sV[s]);
991: } else {
992: const int numSkip = pO-1;
994: sV[s].setSkip(sV[s].getSize()+numSkip);
995: sV[s].setLimit(sV[s].getSize()+pSupportSize);
996: sieve.orientedSupport(point, sV[s]);
997: sieve.orientedSupport(point, sV[s]);
998: sV[s].setSkip(0);
999: sV[s].setLimit(0);
1000: }
1001: }
1002: }
1003: sV[1-s].clear();
1004: }
1005: }
1006: };
1008: namespace IFSieveDef {
1009: template<typename PointType_>
1010: class Sequence {
1011: public:
1012: typedef PointType_ point_type;
1013: class const_iterator {
1014: public:
1015: // Standard iterator typedefs
1016: typedef std::input_iterator_tag iterator_category;
1017: typedef PointType_ value_type;
1018: typedef int difference_type;
1019: typedef value_type* pointer;
1020: typedef value_type& reference;
1021: protected:
1022: const point_type *_data;
1023: int _pos;
1024: public:
1025: const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
1026: virtual ~const_iterator() {};
1027: public:
1028: virtual bool operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
1029: virtual bool operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
1030: virtual const value_type operator*() const {return this->_data[this->_pos];};
1031: virtual const_iterator& operator++() {++this->_pos; return *this;};
1032: virtual const_iterator operator++(int n) {
1033: const_iterator tmp(*this);
1034: ++this->_pos;
1035: return tmp;
1036: };
1037: };
1038: typedef const_iterator iterator;
1039: protected:
1040: const point_type *_data;
1041: int _begin;
1042: int _end;
1043: public:
1044: Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
1045: virtual ~Sequence() {};
1046: public:
1047: virtual iterator begin() const {return iterator(_data, _begin);};
1048: virtual iterator end() const {return iterator(_data, _end);};
1049: size_t size() const {return(_end - _begin);}
1050: void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
1051: };
1052: }
1054: // Interval Final Sieve
1055: // This is just two CSR matrices that give cones and supports
1056: // It is completely static and cannot be resized
1057: // It will operator on visitors, rather than sequences (which are messy)
1058: template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1059: class IFSieve : public ParallelObject {
1060: public:
1061: // Types
1062: typedef IFSieve<Point_,Allocator_> this_type;
1063: typedef Point_ point_type;
1064: typedef SimpleArrow<point_type,point_type> arrow_type;
1065: typedef typename arrow_type::source_type source_type;
1066: typedef typename arrow_type::target_type target_type;
1067: typedef int index_type;
1068: // Allocators
1069: typedef Allocator_ point_allocator_type;
1070: typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1071: typedef typename point_allocator_type::template rebind<int>::other int_allocator_type;
1072: // Interval
1073: typedef Interval<point_type, point_allocator_type> chart_type;
1074: // Dynamic structure
1075: typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1076: // Iterator interface
1077: typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1078: typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1079: // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1080: typedef std::set<point_type> pointSet;
1081: typedef ALE::array<point_type> pointArray;
1082: typedef pointSet coneSet;
1083: typedef pointSet supportSet;
1084: typedef pointArray coneArray;
1085: typedef pointArray supportArray;
1086: protected:
1087: // Arrow Containers
1088: typedef index_type *offsets_type;
1089: typedef point_type *cones_type;
1090: typedef point_type *supports_type;
1091: // Decorators
1092: typedef int *orientations_type;
1093: protected:
1094: // Data
1095: bool indexAllocated;
1096: offsets_type coneOffsets;
1097: offsets_type supportOffsets;
1098: bool pointAllocated;
1099: index_type maxConeSize;
1100: index_type maxSupportSize;
1101: index_type baseSize;
1102: index_type capSize;
1103: cones_type cones;
1104: supports_type supports;
1105: bool orientCones;
1106: orientations_type coneOrientations;
1107: chart_type chart;
1108: int_allocator_type intAlloc;
1109: index_allocator_type indexAlloc;
1110: point_allocator_type pointAlloc;
1111: newpoints_type newCones;
1112: newpoints_type newSupports;
1113: // Sequences
1114: Obj<coneSequence> coneSeq;
1115: Obj<supportSequence> supportSeq;
1116: protected: // Memory Management
1117: void createIndices() {
1118: this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1119: this->coneOffsets -= this->chart.min();
1120: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1121: this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1122: this->supportOffsets -= this->chart.min();
1123: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1124: this->indexAllocated = true;
1125: };
1126: void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1127: if (*coneOffsets) {
1128: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1129: *coneOffsets += chart.min();
1130: indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1131: *coneOffsets = NULL;
1132: }
1133: if (*supportOffsets) {
1134: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1135: *supportOffsets += chart.min();
1136: indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1137: *supportOffsets = NULL;
1138: }
1139: };
1140: void destroyIndices() {
1141: this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1142: this->indexAllocated = false;
1143: this->maxConeSize = -1;
1144: this->maxSupportSize = -1;
1145: this->baseSize = -1;
1146: this->capSize = -1;
1147: };
1148: void createPoints() {
1149: this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1150: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1151: this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1152: for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1153: if (orientCones) {
1154: this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1155: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1156: }
1157: this->pointAllocated = true;
1158: };
1159: void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1160: if (*cones) {
1161: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1162: pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1163: *cones = NULL;
1164: }
1165: if (*supports) {
1166: for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1167: pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1168: *supports = NULL;
1169: }
1170: if (*coneOrientations) {
1171: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1172: intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1173: *coneOrientations = NULL;
1174: }
1175: };
1176: void destroyPoints() {
1177: this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1178: this->pointAllocated = false;
1179: };
1180: void prefixSum(const offsets_type array) {
1181: for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1182: array[p] = array[p] + array[p-1];
1183: }
1184: };
1185: void calculateBaseAndCapSize() {
1186: this->baseSize = 0;
1187: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1188: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1189: ++this->baseSize;
1190: }
1191: }
1192: this->capSize = 0;
1193: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1194: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1195: ++this->capSize;
1196: }
1197: }
1198: };
1199: public:
1200: 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) {
1201: this->coneSeq = new coneSequence(NULL, 0, 0);
1202: this->supportSeq = new supportSequence(NULL, 0, 0);
1203: };
1204: 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) {
1205: this->setChart(chart_type(min, max));
1206: this->coneSeq = new coneSequence(NULL, 0, 0);
1207: this->supportSeq = new supportSequence(NULL, 0, 0);
1208: };
1209: ~IFSieve() {
1210: this->destroyPoints();
1211: this->destroyIndices();
1212: };
1213: public: // Accessors
1214: const chart_type& getChart() const {return this->chart;};
1215: void setChart(const chart_type& chart) {
1216: this->destroyPoints();
1217: this->destroyIndices();
1218: this->chart = chart;
1219: this->createIndices();
1220: };
1221: index_type getMaxConeSize() const {
1222: return this->maxConeSize;
1223: };
1224: index_type getMaxSupportSize() const {
1225: return this->maxSupportSize;
1226: };
1227: bool baseContains(const point_type& p) const {
1228: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1229: this->chart.checkPoint(p);
1231: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1232: return true;
1233: }
1234: return false;
1235: };
1236: bool orientedCones() const {return this->orientCones;};
1237: public: // Construction
1238: index_type getConeSize(const point_type& p) const {
1239: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1240: this->chart.checkPoint(p);
1241: return this->coneOffsets[p+1]-this->coneOffsets[p];
1242: };
1243: void setConeSize(const point_type& p, const index_type c) {
1244: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1245: this->chart.checkPoint(p);
1246: this->coneOffsets[p+1] = c;
1247: this->maxConeSize = std::max(this->maxConeSize, c);
1248: };
1249: void addConeSize(const point_type& p, const index_type c) {
1250: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1251: this->chart.checkPoint(p);
1252: this->coneOffsets[p+1] += c;
1253: this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1254: };
1255: index_type getSupportSize(const point_type& p) const {
1256: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1257: this->chart.checkPoint(p);
1258: return this->supportOffsets[p+1]-this->supportOffsets[p];
1259: };
1260: void setSupportSize(const point_type& p, const index_type s) {
1261: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1262: this->chart.checkPoint(p);
1263: this->supportOffsets[p+1] = s;
1264: this->maxSupportSize = std::max(this->maxSupportSize, s);
1265: };
1266: void addSupportSize(const point_type& p, const index_type s) {
1267: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1268: this->chart.checkPoint(p);
1269: this->supportOffsets[p+1] += s;
1270: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1271: };
1272: void allocate() {
1273: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1274: this->prefixSum(this->coneOffsets);
1275: this->prefixSum(this->supportOffsets);
1276: this->createPoints();
1277: this->calculateBaseAndCapSize();
1278: }
1279: // Purely for backwards compatibility
1280: template<typename Color>
1281: void addArrow(const point_type& p, const point_type& q, const Color c, const bool forceSupport = false) {
1282: this->addArrow(p, q, forceSupport);
1283: }
1284: void addArrow(const point_type& p, const point_type& q, const bool forceSupport = false) {
1285: if (!this->chart.hasPoint(q)) {
1286: if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1287: const index_type start = this->coneOffsets[q];
1288: const index_type end = this->coneOffsets[q+1];
1290: for(int c = start; c < end; ++c) {
1291: this->newCones[q].push_back(this->cones[c]);
1292: }
1293: }
1294: this->newCones[q].push_back(p);
1295: }
1296: if (!this->chart.hasPoint(p) || forceSupport) {
1297: if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1298: const index_type start = this->supportOffsets[p];
1299: const index_type end = this->supportOffsets[p+1];
1301: for(int s = start; s < end; ++s) {
1302: this->newSupports[p].push_back(this->supports[s]);
1303: }
1304: }
1305: this->newSupports[p].push_back(q);
1306: }
1307: };
1308: void reallocate() {
1309: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1310: if (!this->newCones.size() && !this->newSupports.size()) return;
1311: const chart_type oldChart = this->chart;
1312: offsets_type oldConeOffsets = this->coneOffsets;
1313: offsets_type oldSupportOffsets = this->supportOffsets;
1314: cones_type oldCones = this->cones;
1315: supports_type oldSupports = this->supports;
1316: orientations_type oldConeOrientations = this->coneOrientations;
1317: point_type min = this->chart.min();
1318: point_type max = this->chart.max()-1;
1320: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1321: min = std::min(min, c_iter->first);
1322: max = std::max(max, c_iter->first);
1323: }
1324: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1325: min = std::min(min, s_iter->first);
1326: max = std::max(max, s_iter->first);
1327: }
1328: this->chart = chart_type(min, max+1);
1329: this->createIndices();
1330: // Copy sizes (converted from offsets)
1331: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1332: this->coneOffsets[p+1] = oldConeOffsets[p+1]-oldConeOffsets[p];
1333: this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1334: }
1335: // Inject new sizes
1336: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1337: this->coneOffsets[c_iter->first+1] = c_iter->second.size();
1338: this->maxConeSize = std::max(this->maxConeSize, (int) c_iter->second.size());
1339: }
1340: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1341: this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1342: this->maxSupportSize = std::max(this->maxSupportSize, (int) s_iter->second.size());
1343: }
1344: this->prefixSum(this->coneOffsets);
1345: this->prefixSum(this->supportOffsets);
1346: this->createPoints();
1347: this->calculateBaseAndCapSize();
1348: // Copy cones and supports
1349: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1350: const index_type cStart = this->coneOffsets[p];
1351: const index_type cOStart = oldConeOffsets[p];
1352: const index_type cOEnd = oldConeOffsets[p+1];
1353: const index_type sStart = this->supportOffsets[p];
1354: const index_type sOStart = oldSupportOffsets[p];
1355: const index_type sOEnd = oldSupportOffsets[p+1];
1357: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1358: this->cones[c] = oldCones[cO];
1359: }
1360: for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1361: this->supports[s] = oldSupports[sO];
1362: }
1363: if (this->orientCones) {
1364: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1365: this->coneOrientations[c] = oldConeOrientations[cO];
1366: }
1367: }
1368: }
1369: // Inject new cones and supports
1370: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1371: index_type start = this->coneOffsets[c_iter->first];
1373: for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1374: this->cones[start++] = *p_iter;
1375: }
1376: if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1377: }
1378: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1379: index_type start = this->supportOffsets[s_iter->first];
1381: for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1382: this->supports[start++] = *p_iter;
1383: }
1384: if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1385: }
1386: this->newCones.clear();
1387: this->newSupports.clear();
1388: this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1389: this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1390: };
1391: // Recalculate the support offsets and fill the supports
1392: // This is used when an IFSieve is being used as a label
1393: void recalculateLabel() {
1394: ISieveVisitor::PointRetriever<IFSieve> v(1);
1396: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1397: this->supportOffsets[p+1] = 0;
1398: }
1399: this->maxSupportSize = 0;
1400: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1401: this->cone(p, v);
1402: if (v.getSize()) {
1403: this->supportOffsets[v.getPoints()[0]+1]++;
1404: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1405: }
1406: v.clear();
1407: }
1408: this->prefixSum(this->supportOffsets);
1409: this->calculateBaseAndCapSize();
1410: this->symmetrize();
1411: };
1412: void setCone(const point_type cone[], const point_type& p) {
1413: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1414: this->chart.checkPoint(p);
1415: const index_type start = this->coneOffsets[p];
1416: const index_type end = this->coneOffsets[p+1];
1418: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1419: this->cones[c] = cone[i];
1420: }
1421: };
1422: void setCone(const point_type cone, const point_type& p) {
1423: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1424: this->chart.checkPoint(p);
1425: const index_type start = this->coneOffsets[p];
1426: const index_type end = this->coneOffsets[p+1];
1428: if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1429: this->cones[start] = cone;
1430: };
1431: #if 0
1432: template<typename PointSequence>
1433: void setCone(const PointSequence& cone, const point_type& p) {
1434: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1435: this->chart.checkPoint(p);
1436: const index_type start = this->coneOffsets[p];
1437: const index_type end = this->coneOffsets[p+1];
1438: if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1439: typename PointSequence::iterator c_iter = cone.begin();
1441: for(index_type c = start; c < end; ++c, ++c_iter) {
1442: this->cones[c] = *c_iter;
1443: }
1444: };
1445: #endif
1446: void setSupport(const point_type& p, const point_type support[]) {
1447: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1448: this->chart.checkPoint(p);
1449: const index_type start = this->supportOffsets[p];
1450: const index_type end = this->supportOffsets[p+1];
1452: for(index_type s = start, i = 0; s < end; ++s, ++i) {
1453: this->supports[s] = support[i];
1454: }
1455: };
1456: #if 0
1457: template<typename PointSequence>
1458: void setSupport(const point_type& p, const PointSequence& support) {
1459: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1460: this->chart.checkPoint(p);
1461: const index_type start = this->supportOffsets[p];
1462: const index_type end = this->supportOffsets[p+1];
1463: if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1464: typename PointSequence::iterator s_iter = support.begin();
1466: for(index_type s = start; s < end; ++s, ++s_iter) {
1467: this->supports[s] = *s_iter;
1468: }
1469: };
1470: #endif
1471: void setConeOrientation(const int coneOrientation[], const point_type& p) {
1472: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1473: this->chart.checkPoint(p);
1474: const index_type start = this->coneOffsets[p];
1475: const index_type end = this->coneOffsets[p+1];
1477: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1478: this->coneOrientations[c] = coneOrientation[i];
1479: }
1480: };
1481: void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1482: for(point_type p = 0; p < numCells; ++p) {
1483: const index_type start = p*numCorners;
1484: const index_type end = (p+1)*numCorners;
1486: for(index_type c = start; c < end; ++c) {
1487: const point_type q = cones[c]+offset;
1489: this->supportOffsets[q+1]++;
1490: }
1491: }
1492: for(point_type p = numCells; p < this->chart.max(); ++p) {
1493: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1494: }
1495: };
1496: void symmetrize() {
1497: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1498: offsets -= this->chart.min();
1499: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1500: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1502: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1503: const index_type start = this->coneOffsets[p];
1504: const index_type end = this->coneOffsets[p+1];
1506: for(index_type c = start; c < end; ++c) {
1507: const point_type q = this->cones[c];
1509: this->chart.checkPoint(q);
1510: this->supports[this->supportOffsets[q]+offsets[q]] = p;
1511: ++offsets[q];
1512: }
1513: }
1514: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1515: indexAlloc.deallocate(offsets, this->chart.size()+1);
1516: }
1517: index_type getBaseSize() const {
1518: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1519: return this->baseSize;
1520: }
1521: index_type getCapSize() const {
1522: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1523: return this->capSize;
1524: }
1525: public: // Traversals
1526: template<typename Visitor>
1527: void roots(const Visitor& v) const {
1528: this->roots(const_cast<Visitor&>(v));
1529: }
1530: template<typename Visitor>
1531: void roots(Visitor& v) const {
1532: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1534: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1535: if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1536: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1537: v.visitPoint(p);
1538: }
1539: }
1540: }
1541: }
1542: template<typename Visitor>
1543: void leaves(const Visitor& v) const {
1544: this->leaves(const_cast<Visitor&>(v));
1545: }
1546: template<typename Visitor>
1547: void leaves(Visitor& v) const {
1548: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1550: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1551: if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1552: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1553: v.visitPoint(p);
1554: }
1555: }
1556: }
1557: }
1558: template<typename Visitor>
1559: void base(const Visitor& v) const {
1560: this->base(const_cast<Visitor&>(v));
1561: }
1562: template<typename Visitor>
1563: void base(Visitor& v) const {
1564: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1566: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1567: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1568: v.visitPoint(p);
1569: }
1570: }
1571: }
1572: template<typename Visitor>
1573: void cap(const Visitor& v) const {
1574: this->cap(const_cast<Visitor&>(v));
1575: }
1576: template<typename Visitor>
1577: void cap(Visitor& v) const {
1578: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1580: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1581: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1582: v.visitPoint(p);
1583: }
1584: }
1585: }
1586: template<typename PointSequence, typename Visitor>
1587: void cone(const PointSequence& points, Visitor& v) const {
1588: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1589: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1590: const point_type p = *p_iter;
1591: this->chart.checkPoint(p);
1592: const index_type start = this->coneOffsets[p];
1593: const index_type end = this->coneOffsets[p+1];
1595: for(index_type c = start; c < end; ++c) {
1596: v.visitArrow(arrow_type(this->cones[c], p));
1597: v.visitPoint(this->cones[c]);
1598: }
1599: }
1600: }
1601: template<typename Visitor>
1602: void cone(const point_type& p, Visitor& v) const {
1603: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1604: this->chart.checkPoint(p);
1605: const index_type start = this->coneOffsets[p];
1606: const index_type end = this->coneOffsets[p+1];
1608: for(index_type c = start; c < end; ++c) {
1609: v.visitArrow(arrow_type(this->cones[c], p));
1610: v.visitPoint(this->cones[c]);
1611: }
1612: }
1613: const Obj<coneSequence>& cone(const point_type& p) const {
1614: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1615: if (!this->chart.hasPoint(p)) {
1616: this->coneSeq->reset(this->cones, 0, 0);
1617: } else {
1618: this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1619: }
1620: return this->coneSeq;
1621: }
1622: template<typename PointSequence, typename Visitor>
1623: void support(const PointSequence& points, Visitor& v) const {
1624: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1625: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1626: const point_type p = *p_iter;
1627: this->chart.checkPoint(p);
1628: const index_type start = this->supportOffsets[p];
1629: const index_type end = this->supportOffsets[p+1];
1631: for(index_type s = start; s < end; ++s) {
1632: v.visitArrow(arrow_type(p, this->supports[s]));
1633: v.visitPoint(this->supports[s]);
1634: }
1635: }
1636: }
1637: template<typename Visitor>
1638: void support(const point_type& p, Visitor& v) const {
1639: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1640: this->chart.checkPoint(p);
1641: const index_type start = this->supportOffsets[p];
1642: const index_type end = this->supportOffsets[p+1];
1644: for(index_type s = start; s < end; ++s) {
1645: v.visitArrow(arrow_type(p, this->supports[s]));
1646: v.visitPoint(this->supports[s]);
1647: }
1648: }
1649: const Obj<supportSequence>& support(const point_type& p) const {
1650: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1651: if (!this->chart.hasPoint(p)) {
1652: this->supportSeq->reset(this->supports, 0, 0);
1653: } else {
1654: this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1655: }
1656: return this->supportSeq;
1657: }
1658: template<typename Visitor>
1659: void orientedCone(const point_type& p, Visitor& v) const {
1660: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1661: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1662: this->chart.checkPoint(p);
1663: const index_type start = this->coneOffsets[p];
1664: const index_type end = this->coneOffsets[p+1];
1666: for(index_type c = start; c < end; ++c) {
1667: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1668: v.visitPoint(this->cones[c], this->coneOrientations[c]);
1669: }
1670: }
1671: template<typename Visitor>
1672: void orientedReverseCone(const point_type& p, Visitor& v) const {
1673: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1674: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1675: this->chart.checkPoint(p);
1676: const index_type start = this->coneOffsets[p];
1677: const index_type end = this->coneOffsets[p+1];
1679: for(index_type c = end-1; c >= start; --c) {
1680: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1681: v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1682: }
1683: }
1684: template<typename Visitor>
1685: void orientedSupport(const point_type& p, Visitor& v) const {
1686: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1687: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1688: this->chart.checkPoint(p);
1689: const index_type start = this->supportOffsets[p];
1690: const index_type end = this->supportOffsets[p+1];
1692: for(index_type s = start; s < end; ++s) {
1693: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1694: //v.visitPoint(this->supports[s], this->supportOrientations[s]);
1695: v.visitArrow(arrow_type(this->supports[s], p), 0);
1696: v.visitPoint(this->supports[s], 0);
1697: }
1698: }
1699: template<typename Visitor>
1700: void orientedReverseSupport(const point_type& p, Visitor& v) const {
1701: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1702: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1703: this->chart.checkPoint(p);
1704: const index_type start = this->supportOffsets[p];
1705: const index_type end = this->supportOffsets[p+1];
1707: for(index_type s = end-1; s >= start; --s) {
1708: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1709: //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
1710: v.visitArrow(arrow_type(this->supports[s], p), 0);
1711: v.visitPoint(this->supports[s], 0);
1712: }
1713: }
1714: // Currently does only 1 level
1715: // Does not check for uniqueness
1716: template<typename Visitor>
1717: void meet(const point_type& p, const point_type& q, Visitor& v) const {
1718: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1719: this->chart.checkPoint(p);
1720: this->chart.checkPoint(q);
1721: const index_type startP = this->coneOffsets[p];
1722: const index_type endP = this->coneOffsets[p+1];
1723: const index_type startQ = this->coneOffsets[q];
1724: const index_type endQ = this->coneOffsets[q+1];
1726: for(index_type cP = startP; cP < endP; ++cP) {
1727: const point_type& c1 = this->cones[cP];
1729: for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1730: if (c1 == this->cones[cQ]) v.visitPoint(c1);
1731: }
1732: if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1733: }
1734: }
1735: // Currently does only 1 level
1736: template<typename Sequence, typename Visitor>
1737: void join(const Sequence& points, Visitor& v) const {
1738: typedef std::set<point_type> pointSet;
1739: pointSet intersect[2] = {pointSet(), pointSet()};
1740: pointSet tmp;
1741: int p = 0;
1742: int c = 0;
1744: for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1745: this->chart.checkPoint(*p_iter);
1746: tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1747: if (p == 0) {
1748: intersect[1-c].insert(tmp.begin(), tmp.end());
1749: p++;
1750: } else {
1751: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1752: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1753: intersect[c].clear();
1754: }
1755: c = 1 - c;
1756: tmp.clear();
1757: }
1758: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1759: v.visitPoint(*p_iter);
1760: }
1761: }
1762: public: // Viewing
1763: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1764: ostringstream txt;
1765: int rank;
1767: if (comm == MPI_COMM_NULL) {
1768: comm = this->comm();
1769: rank = this->commRank();
1770: } else {
1771: MPI_Comm_rank(comm, &rank);
1772: }
1773: if (name == "") {
1774: if(rank == 0) {
1775: txt << "viewing an IFSieve" << std::endl;
1776: }
1777: } else {
1778: if(rank == 0) {
1779: txt << "viewing IFSieve '" << name << "'" << std::endl;
1780: }
1781: }
1782: if(rank == 0) {
1783: txt << "cap --> base:" << std::endl;
1784: }
1785: ISieveVisitor::PrintVisitor pV(txt, rank);
1786: this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1787: PetscSynchronizedPrintf(comm, txt.str().c_str());
1788: PetscSynchronizedFlush(comm);
1789: ostringstream txt2;
1791: if(rank == 0) {
1792: txt2 << "base <-- cap:" << std::endl;
1793: }
1794: ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1795: this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1796: PetscSynchronizedPrintf(comm, txt2.str().c_str());
1797: PetscSynchronizedFlush(comm);
1798: if (orientCones) {
1799: ostringstream txt3;
1801: if(rank == 0) {
1802: txt3 << "Orientation:" << std::endl;
1803: }
1804: ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1805: this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1806: PetscSynchronizedPrintf(comm, txt3.str().c_str());
1807: PetscSynchronizedFlush(comm);
1808: }
1809: };
1810: };
1812: class ISieveConverter {
1813: public:
1814: template<typename Sieve, typename ISieve, typename Renumbering>
1815: static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1816: // First construct a renumbering of the sieve points
1817: const Obj<typename Sieve::baseSequence>& base = sieve.base();
1818: const Obj<typename Sieve::capSequence>& cap = sieve.cap();
1819: typename ISieve::point_type min = 0;
1820: typename ISieve::point_type max = 0;
1822: if (renumber) {
1823: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1824: renumbering[*b_iter] = max++;
1825: }
1826: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1827: if (renumbering.find(*c_iter) == renumbering.end()) {
1828: renumbering[*c_iter] = max++;
1829: }
1830: }
1831: } else {
1832: if (base->size()) {
1833: min = *base->begin();
1834: max = *base->begin();
1835: } else if (cap->size()) {
1836: min = *cap->begin();
1837: max = *cap->begin();
1838: }
1839: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1840: min = std::min(min, *b_iter);
1841: max = std::max(max, *b_iter);
1842: }
1843: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1844: min = std::min(min, *c_iter);
1845: max = std::max(max, *c_iter);
1846: }
1847: if (base->size() || cap->size()) {
1848: ++max;
1849: }
1850: for(typename ISieve::point_type p = min; p < max; ++p) {
1851: renumbering[p] = p;
1852: }
1853: }
1854: // Create the ISieve
1855: isieve.setChart(typename ISieve::chart_type(min, max));
1856: // Set cone and support sizes
1857: size_t maxSize = 0;
1859: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1860: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1862: isieve.setConeSize(renumbering[*b_iter], cone->size());
1863: maxSize = std::max(maxSize, cone->size());
1864: }
1865: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1866: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1868: isieve.setSupportSize(renumbering[*c_iter], support->size());
1869: maxSize = std::max(maxSize, support->size());
1870: }
1871: isieve.allocate();
1872: // Fill up cones and supports
1873: typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];
1875: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1876: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1877: int i = 0;
1879: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1880: points[i] = renumbering[*c_iter];
1881: }
1882: isieve.setCone(points, renumbering[*b_iter]);
1883: }
1884: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1885: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
1886: int i = 0;
1888: for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
1889: points[i] = renumbering[*s_iter];
1890: }
1891: isieve.setSupport(renumbering[*c_iter], points);
1892: }
1893: delete [] points;
1894: }
1895: template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
1896: static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
1897: if (isieve.getMaxConeSize() < 0) return;
1898: const Obj<typename Sieve::baseSequence>& base = sieve.base();
1899: int *orientations = new int[isieve.getMaxConeSize()];
1901: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1902: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
1903: int i = 0;
1905: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
1906: typename ArrowSection::point_type arrow(*c_iter, *b_iter);
1908: orientations[i] = orientation->restrictPoint(arrow)[0];
1909: }
1910: isieve.setConeOrientation(orientations, renumbering[*b_iter]);
1911: }
1912: delete [] orientations;
1913: }
1914: template<typename Section, typename ISection, typename Renumbering>
1915: static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
1916: const typename Section::chart_type& chart = coordinates.getChart();
1917: typename ISection::point_type min = *chart.begin();
1918: typename ISection::point_type max = *chart.begin();
1920: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1921: min = std::min(min, *p_iter);
1922: max = std::max(max, *p_iter);
1923: }
1924: icoordinates.setChart(typename ISection::chart_type(min, max+1));
1925: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1926: icoordinates.setFiberDimension(*p_iter, coordinates.getFiberDimension(*p_iter));
1927: }
1928: icoordinates.allocatePoint();
1929: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
1930: icoordinates.updatePoint(*p_iter, coordinates.restrictPoint(*p_iter));
1931: }
1932: }
1933: template<typename IMesh, typename Label>
1934: static void convertLabel(IMesh& imesh, const std::string& name, const Obj<Label>& oldLabel) {
1935: const Obj<typename IMesh::label_type>& label = imesh.createLabel(name);
1936: const typename IMesh::sieve_type::chart_type& chart = imesh.getSieve()->getChart();
1937: int size = 0;
1939: label->setChart(chart);
1940: for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1941: const int coneSize = oldLabel->cone(p)->size();
1943: label->setConeSize(p, coneSize);
1944: size += coneSize;
1945: }
1946: if (size) {label->setSupportSize(0, size);}
1947: label->allocate();
1948: for(typename IMesh::point_type p = chart.min(); p < chart.max(); ++p) {
1949: const Obj<typename Label::coneSequence>& cone = oldLabel->cone(p);
1951: if (cone->size()) {
1952: label->setCone(*cone->begin(), p);
1953: }
1954: }
1955: label->recalculateLabel();
1956: }
1957: template<typename Mesh, typename IMesh, typename Renumbering>
1958: static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
1959: convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
1960: imesh.stratify();
1961: convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
1962: convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
1963: if (mesh.hasRealSection("normals")) {
1964: convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
1965: }
1966: const typename Mesh::labels_type& labels = mesh.getLabels();
1967: std::string heightName("height");
1968: std::string depthName("depth");
1970: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
1971: #ifdef IMESH_NEW_LABELS
1972: if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
1973: convertLabel(imesh, l_iter->first, l_iter->second);
1974: }
1975: #else
1976: imesh.setLabel(l_iter->first, l_iter->second);
1977: #endif
1978: }
1979: }
1980: };
1982: class ISieveSerializer {
1983: public:
1984: template<typename ISieve>
1985: static void writeSieve(const std::string& filename, ISieve& sieve) {
1986: std::ofstream fs;
1988: if (sieve.commRank() == 0) {
1989: fs.open(filename.c_str());
1990: }
1991: writeSieve(fs, sieve);
1992: if (sieve.commRank() == 0) {
1993: fs.close();
1994: }
1995: };
1996: template<typename ISieve>
1997: static void writeSieve(std::ofstream& fs, ISieve& sieve) {
1998: typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
1999: const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
2000: typename ISieve::point_type min = chart->min();
2001: typename ISieve::point_type max = chart->max();
2002: PetscInt *mins = new PetscInt[sieve.commSize()];
2003: PetscInt *maxs = new PetscInt[sieve.commSize()];
2005: // Write sizes
2006: if (sieve.commRank() == 0) {
2007: // Write local
2008: fs << min <<" "<< max << std::endl;
2009: for(typename ISieve::point_type p = min; p < max; ++p) {
2010: fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
2011: }
2012: // Receive and write remote
2013: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2014: PetscInt minp, maxp;
2015: PetscInt *sizes;
2016: MPI_Status status;
2019: MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2020: MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2021: PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2022: MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2023: fs << minp <<" "<< maxp << std::endl;
2024: for(PetscInt s = 0; s < maxp - minp; ++s) {
2025: fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
2026: }
2027: PetscFree(sizes);CHKERRXX(ierr);
2028: mins[pr] = minp;
2029: maxs[pr] = maxp;
2030: }
2031: } else {
2032: // Send remote
2033: //PetscInt min = chart->min();
2034: //PetscInt max = chart->max();
2035: PetscInt s = 0;
2036: PetscInt *sizes;
2039: MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2040: MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2041: PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
2042: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2043: sizes[s*2+0] = sieve.getConeSize(p);
2044: sizes[s*2+1] = sieve.getSupportSize(p);
2045: }
2046: MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2047: PetscFree(sizes);CHKERRXX(ierr);
2048: }
2049: // Write covering
2050: if (sieve.commRank() == 0) {
2051: // Write local
2052: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2054: for(typename ISieve::point_type p = min; p < max; ++p) {
2055: sieve.cone(p, pV);
2056: const typename Visitor::point_type *cone = pV.getPoints();
2057: const int cSize = pV.getSize();
2059: if (cSize > 0) {
2060: for(int c = 0; c < cSize; ++c) {
2061: if (c) {fs << " ";}
2062: fs << cone[c];
2063: }
2064: fs << std::endl;
2065: }
2066: pV.clear();
2068: sieve.orientedCone(p, pV);
2069: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2070: const int oSize = pV.getOrientedSize();
2072: if (oSize > 0) {
2073: for(int c = 0; c < oSize; ++c) {
2074: if (c) {fs << " ";}
2075: fs << oCone[c].second;
2076: }
2077: fs << std::endl;
2078: }
2079: pV.clear();
2081: sieve.support(p, pV);
2082: const typename Visitor::point_type *support = pV.getPoints();
2083: const int sSize = pV.getSize();
2085: if (sSize > 0) {
2086: for(int s = 0; s < sSize; ++s) {
2087: if (s) {fs << " ";}
2088: fs << support[s];
2089: }
2090: fs << std::endl;
2091: }
2092: pV.clear();
2093: }
2094: // Receive and write remote
2095: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2096: PetscInt size;
2097: PetscInt *data;
2098: PetscInt off = 0;
2099: MPI_Status status;
2102: MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2103: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2104: MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2105: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2106: PetscInt cSize = data[off++];
2108: fs << cSize << std::endl;
2109: if (cSize > 0) {
2110: for(int c = 0; c < cSize; ++c) {
2111: if (c) {fs << " ";}
2112: fs << data[off++];
2113: }
2114: fs << std::endl;
2115: }
2116: PetscInt oSize = data[off++];
2118: if (oSize > 0) {
2119: for(int c = 0; c < oSize; ++c) {
2120: if (c) {fs << " ";}
2121: fs << data[off++];
2122: }
2123: fs << std::endl;
2124: }
2125: PetscInt sSize = data[off++];
2127: fs << sSize << std::endl;
2128: if (sSize > 0) {
2129: for(int s = 0; s < sSize; ++s) {
2130: if (s) {fs << " ";}
2131: fs << data[off++];
2132: }
2133: fs << std::endl;
2134: }
2135: }
2136: assert(off == size);
2137: PetscFree(data);CHKERRXX(ierr);
2138: }
2139: } else {
2140: // Send remote
2141: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2142: PetscInt totalConeSize = 0;
2143: PetscInt totalSupportSize = 0;
2145: for(typename ISieve::point_type p = min; p < max; ++p) {
2146: totalConeSize += sieve.getConeSize(p);
2147: totalSupportSize += sieve.getSupportSize(p);
2148: }
2149: PetscInt size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
2150: PetscInt off = 0;
2151: PetscInt *data;
2154: MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2155: // There is no nice way to make a generic MPI type here. Really sucky
2156: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2157: for(typename ISieve::point_type p = min; p < max; ++p) {
2158: sieve.cone(p, pV);
2159: const typename Visitor::point_type *cone = pV.getPoints();
2160: const int cSize = pV.getSize();
2162: data[off++] = cSize;
2163: for(int c = 0; c < cSize; ++c) {
2164: data[off++] = cone[c];
2165: }
2166: pV.clear();
2168: sieve.orientedCone(p, pV);
2169: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2170: const int oSize = pV.getOrientedSize();
2172: data[off++] = oSize;
2173: for(int c = 0; c < oSize; ++c) {
2174: data[off++] = oCone[c].second;
2175: }
2176: pV.clear();
2178: sieve.support(p, pV);
2179: const typename Visitor::point_type *support = pV.getPoints();
2180: const int sSize = pV.getSize();
2182: data[off++] = sSize;
2183: for(int s = 0; s < sSize; ++s) {
2184: data[off++] = support[s];
2185: }
2186: pV.clear();
2187: }
2188: MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2189: PetscFree(data);CHKERRXX(ierr);
2190: }
2191: delete [] mins;
2192: delete [] maxs;
2193: // Output renumbering
2194: };
2195: template<typename ISieve>
2196: static void loadSieve(const std::string& filename, ISieve& sieve) {
2197: std::ifstream fs;
2199: if (sieve.commRank() == 0) {
2200: fs.open(filename.c_str());
2201: }
2202: loadSieve(fs, sieve);
2203: if (sieve.commRank() == 0) {
2204: fs.close();
2205: }
2206: };
2207: template<typename ISieve>
2208: static void loadSieve(std::ifstream& fs, ISieve& sieve) {
2209: typename ISieve::point_type min, max;
2210: PetscInt *mins = new PetscInt[sieve.commSize()];
2211: PetscInt *maxs = new PetscInt[sieve.commSize()];
2212: PetscInt *totalConeSizes = new PetscInt[sieve.commSize()];
2213: PetscInt *totalSupportSizes = new PetscInt[sieve.commSize()];
2215: // Load sizes
2216: if (sieve.commRank() == 0) {
2217: // Load local
2218: fs >> min;
2219: fs >> max;
2220: sieve.setChart(typename ISieve::chart_type(min, max));
2221: for(typename ISieve::point_type p = min; p < max; ++p) {
2222: typename ISieve::index_type coneSize, supportSize;
2224: fs >> coneSize;
2225: fs >> supportSize;
2226: sieve.setConeSize(p, coneSize);
2227: sieve.setSupportSize(p, supportSize);
2228: }
2229: // Load and send remote
2230: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2231: PetscInt minp, maxp;
2232: PetscInt *sizes;
2235: fs >> minp;
2236: fs >> maxp;
2237: MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2238: MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2239: PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2240: totalConeSizes[pr] = 0;
2241: totalSupportSizes[pr] = 0;
2242: for(PetscInt s = 0; s < maxp - minp; ++s) {
2243: fs >> sizes[s*2+0];
2244: fs >> sizes[s*2+1];
2245: totalConeSizes[pr] += sizes[s*2+0];
2246: totalSupportSizes[pr] += sizes[s*2+1];
2247: }
2248: MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2249: PetscFree(sizes);CHKERRXX(ierr);
2250: mins[pr] = minp;
2251: maxs[pr] = maxp;
2252: }
2253: } else {
2254: // Load remote
2255: PetscInt s = 0;
2256: PetscInt *sizes;
2257: MPI_Status status;
2260: MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2261: MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2262: sieve.setChart(typename ISieve::chart_type(min, max));
2263: PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2264: MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2265: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2266: sieve.setConeSize(p, sizes[s*2+0]);
2267: sieve.setSupportSize(p, sizes[s*2+1]);
2268: }
2269: PetscFree(sizes);CHKERRXX(ierr);
2270: }
2271: sieve.allocate();
2272: // Load covering
2273: if (sieve.commRank() == 0) {
2274: // Load local
2275: typename ISieve::index_type maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
2276: typename ISieve::point_type *points = new typename ISieve::point_type[maxSize];
2278: for(typename ISieve::point_type p = min; p < max; ++p) {
2279: typename ISieve::index_type coneSize = sieve.getConeSize(p);
2280: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2282: if (coneSize > 0) {
2283: for(int c = 0; c < coneSize; ++c) {
2284: fs >> points[c];
2285: }
2286: sieve.setCone(points, p);
2287: if (sieve.orientedCones()) {
2288: for(int c = 0; c < coneSize; ++c) {
2289: fs >> points[c];
2290: }
2291: sieve.setConeOrientation(points, p);
2292: }
2293: }
2294: if (supportSize > 0) {
2295: for(int s = 0; s < supportSize; ++s) {
2296: fs >> points[s];
2297: }
2298: sieve.setSupport(p, points);
2299: }
2300: }
2301: delete [] points;
2302: // Load and send remote
2303: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2304: PetscInt size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
2305: PetscInt off = 0;
2306: PetscInt *data;
2309: MPI_Send(&size, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2310: // There is no nice way to make a generic MPI type here. Really sucky
2311: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2312: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2313: PetscInt coneSize, supportSize;
2315: fs >> coneSize;
2316: data[off++] = coneSize;
2317: if (coneSize > 0) {
2318: for(int c = 0; c < coneSize; ++c) {
2319: fs >> data[off++];
2320: }
2321: for(int c = 0; c < coneSize; ++c) {
2322: fs >> data[off++];
2323: }
2324: }
2325: fs >> supportSize;
2326: data[off++] = supportSize;
2327: if (supportSize > 0) {
2328: for(int s = 0; s < supportSize; ++s) {
2329: fs >> data[off++];
2330: }
2331: }
2332: }
2333: assert(off == size);
2334: MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2335: PetscFree(data);CHKERRXX(ierr);
2336: }
2337: delete [] mins;
2338: delete [] maxs;
2339: } else {
2340: // Load remote
2341: PetscInt size;
2342: PetscInt *data;
2343: PetscInt off = 0;
2344: MPI_Status status;
2347: MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2348: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2349: MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2350: for(typename ISieve::point_type p = min; p < max; ++p) {
2351: typename ISieve::index_type coneSize = sieve.getConeSize(p);
2352: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2353: PetscInt cs = data[off++];
2355: assert(cs == coneSize);
2356: if (coneSize > 0) {
2357: sieve.setCone(&data[off], p);
2358: off += coneSize;
2359: if (sieve.orientedCones()) {
2360: sieve.setConeOrientation(&data[off], p);
2361: off += coneSize;
2362: }
2363: }
2364: PetscInt ss = data[off++];
2366: assert(ss == supportSize);
2367: if (supportSize > 0) {
2368: sieve.setSupport(p, &data[off]);
2369: off += supportSize;
2370: }
2371: }
2372: assert(off == size);
2373: }
2374: // Load renumbering
2375: };
2376: };
2377: }
2379: #endif