Actual source code: ISieve.hh
1: #ifndef included_ALE_ISieve_hh
2: #define included_ALE_ISieve_hh
4: #ifndef included_ALE_hh
5: #include <sieve/ALE.hh>
6: #endif
8: #include <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: typedef ISieveVisitor::PointRetriever<Sieve,Retriever> TmpRetriever;
801: Retriever pV(200, v, true); // Correct estimate is pow(std::max(1, sieve->getMaxConeSize()), mesh->depth())
802: TmpRetriever cV[2] = {TmpRetriever(200,pV), TmpRetriever(200,pV)};
803: int c = 0;
805: v.visitPoint(p, 0);
806: // Cone is guarateed to be ordered correctly
807: sieve.orientedCone(p, cV[c]);
809: while(cV[c].getOrientedSize()) {
810: const typename Retriever::oriented_point_type *cone = cV[c].getOrientedPoints();
811: const int coneSize = cV[c].getOrientedSize();
812: c = 1 - c;
814: for(int p = 0; p < coneSize; ++p) {
815: const typename Retriever::point_type& point = cone[p].first;
816: int pO = cone[p].second == 0 ? 1 : cone[p].second;
817: const int pConeSize = sieve.getConeSize(point);
819: if (pO < 0) {
820: if (pO == -pConeSize) {
821: sieve.orientedReverseCone(point, cV[c]);
822: } else {
823: const int numSkip = sieve.getConeSize(point) + pO;
825: cV[c].setSkip(cV[c].getSize()+numSkip);
826: cV[c].setLimit(cV[c].getSize()+pConeSize);
827: sieve.orientedReverseCone(point, cV[c]);
828: sieve.orientedReverseCone(point, cV[c]);
829: cV[c].setSkip(0);
830: cV[c].setLimit(0);
831: }
832: } else {
833: if (pO == 1) {
834: sieve.orientedCone(point, cV[c]);
835: } else {
836: const int numSkip = pO-1;
838: cV[c].setSkip(cV[c].getSize()+numSkip);
839: cV[c].setLimit(cV[c].getSize()+pConeSize);
840: sieve.orientedCone(point, cV[c]);
841: sieve.orientedCone(point, cV[c]);
842: cV[c].setSkip(0);
843: cV[c].setLimit(0);
844: }
845: }
846: }
847: cV[1-c].clear();
848: }
849: #if 0
850: // These contain arrows paired with orientations from the \emph{previous} arrow
851: Obj<orientedArrowArray> cone = new orientedArrowArray();
852: Obj<orientedArrowArray> base = new orientedArrowArray();
853: coneSet seen;
855: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
856: const arrow_type arrow(*c_iter, p);
858: cone->push_back(oriented_arrow_type(arrow, 1)); // Notice the orientation of leaf cells is always 1
859: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0])); // However, we return the actual arrow orientation
860: }
861: for(int i = 1; i < depth; ++i) {
862: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
864: cone->clear();
865: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
866: const arrow_type& arrow = b_iter->first; // This is an arrow considered in the previous round
867: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source); // We are going to get the cone of this guy
868: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0]; // The orientation of arrow, which is our pO
870: if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
871: o = -(o+1);
872: }
873: if (o < 0) {
874: const int size = pCone->size();
876: if (o == -size) {
877: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
878: if (seen.find(*c_iter) == seen.end()) {
879: const arrow_type newArrow(*c_iter, arrow.source);
880: int pointO = orientation->restrictPoint(newArrow)[0];
882: seen.insert(*c_iter);
883: cone->push_back(oriented_arrow_type(newArrow, o));
884: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
885: }
886: }
887: } else {
888: const int numSkip = size + o;
889: int count = 0;
891: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
892: if (count < numSkip) continue;
893: if (seen.find(*c_iter) == seen.end()) {
894: const arrow_type newArrow(*c_iter, arrow.source);
895: int pointO = orientation->restrictPoint(newArrow)[0];
897: seen.insert(*c_iter);
898: cone->push_back(oriented_arrow_type(newArrow, o));
899: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
900: }
901: }
902: count = 0;
903: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
904: if (count >= numSkip) break;
905: if (seen.find(*c_iter) == seen.end()) {
906: const arrow_type newArrow(*c_iter, arrow.source);
907: int pointO = orientation->restrictPoint(newArrow)[0];
909: seen.insert(*c_iter);
910: cone->push_back(oriented_arrow_type(newArrow, o));
911: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
912: }
913: }
914: }
915: } else {
916: if (o == 1) {
917: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
918: if (seen.find(*c_iter) == seen.end()) {
919: const arrow_type newArrow(*c_iter, arrow.source);
921: seen.insert(*c_iter);
922: cone->push_back(oriented_arrow_type(newArrow, o));
923: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
924: }
925: }
926: } else {
927: const int numSkip = o-1;
928: int count = 0;
930: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
931: if (count < numSkip) continue;
932: if (seen.find(*c_iter) == seen.end()) {
933: const arrow_type newArrow(*c_iter, arrow.source);
935: seen.insert(*c_iter);
936: cone->push_back(oriented_arrow_type(newArrow, o));
937: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
938: }
939: }
940: count = 0;
941: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
942: if (count >= numSkip) break;
943: if (seen.find(*c_iter) == seen.end()) {
944: const arrow_type newArrow(*c_iter, arrow.source);
946: seen.insert(*c_iter);
947: cone->push_back(oriented_arrow_type(newArrow, o));
948: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
949: }
950: }
951: }
952: }
953: }
954: }
955: #endif
956: }
957: template<typename Visitor>
958: static void orientedStar(const Sieve& sieve, const point_type& p, Visitor& v) {
959: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
960: Retriever sV[2] = {Retriever(200,v), Retriever(200,v)};
961: int s = 0;
963: v.visitPoint(p, 0);
964: // Support is guarateed to be ordered correctly
965: sieve.orientedSupport(p, sV[s]);
967: while(sV[s].getOrientedSize()) {
968: const typename Retriever::oriented_point_type *support = sV[s].getOrientedPoints();
969: const int supportSize = sV[s].getOrientedSize();
970: s = 1 - s;
972: for(int p = 0; p < supportSize; ++p) {
973: const typename Retriever::point_type& point = support[p].first;
974: int pO = support[p].second == 0 ? 1 : support[p].second;
975: const int pSupportSize = sieve.getSupportSize(point);
977: if (pO < 0) {
978: if (pO == -pSupportSize) {
979: sieve.orientedReverseSupport(point, sV[s]);
980: } else {
981: const int numSkip = sieve.getSupportSize(point) + pO;
983: sV[s].setSkip(sV[s].getSize()+numSkip);
984: sV[s].setLimit(sV[s].getSize()+pSupportSize);
985: sieve.orientedReverseSupport(point, sV[s]);
986: sieve.orientedReverseSupport(point, sV[s]);
987: sV[s].setSkip(0);
988: sV[s].setLimit(0);
989: }
990: } else {
991: if (pO == 1) {
992: sieve.orientedSupport(point, sV[s]);
993: } else {
994: const int numSkip = pO-1;
996: sV[s].setSkip(sV[s].getSize()+numSkip);
997: sV[s].setLimit(sV[s].getSize()+pSupportSize);
998: sieve.orientedSupport(point, sV[s]);
999: sieve.orientedSupport(point, sV[s]);
1000: sV[s].setSkip(0);
1001: sV[s].setLimit(0);
1002: }
1003: }
1004: }
1005: sV[1-s].clear();
1006: }
1007: }
1008: };
1010: namespace IFSieveDef {
1011: template<typename PointType_>
1012: class Sequence {
1013: public:
1014: typedef PointType_ point_type;
1015: class const_iterator {
1016: public:
1017: // Standard iterator typedefs
1018: typedef std::input_iterator_tag iterator_category;
1019: typedef PointType_ value_type;
1020: typedef int difference_type;
1021: typedef value_type* pointer;
1022: typedef value_type& reference;
1023: protected:
1024: const point_type *_data;
1025: int _pos;
1026: public:
1027: const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
1028: virtual ~const_iterator() {};
1029: public:
1030: virtual bool operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
1031: virtual bool operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
1032: virtual const value_type operator*() const {return this->_data[this->_pos];};
1033: virtual const_iterator& operator++() {++this->_pos; return *this;};
1034: virtual const_iterator operator++(int n) {
1035: const_iterator tmp(*this);
1036: ++this->_pos;
1037: return tmp;
1038: };
1039: };
1040: typedef const_iterator iterator;
1041: protected:
1042: const point_type *_data;
1043: int _begin;
1044: int _end;
1045: public:
1046: Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
1047: virtual ~Sequence() {};
1048: public:
1049: virtual iterator begin() const {return iterator(_data, _begin);};
1050: virtual iterator end() const {return iterator(_data, _end);};
1051: size_t size() const {return(_end - _begin);}
1052: void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
1053: };
1054: }
1056: // Interval Final Sieve
1057: // This is just two CSR matrices that give cones and supports
1058: // It is completely static and cannot be resized
1059: // It will operator on visitors, rather than sequences (which are messy)
1060: template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1061: class IFSieve : public ParallelObject {
1062: public:
1063: // Types
1064: typedef IFSieve<Point_,Allocator_> this_type;
1065: typedef Point_ point_type;
1066: typedef SimpleArrow<point_type,point_type> arrow_type;
1067: typedef typename arrow_type::source_type source_type;
1068: typedef typename arrow_type::target_type target_type;
1069: typedef int index_type;
1070: // Allocators
1071: typedef Allocator_ point_allocator_type;
1072: typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1073: typedef typename point_allocator_type::template rebind<int>::other int_allocator_type;
1074: // Interval
1075: typedef Interval<point_type, point_allocator_type> chart_type;
1076: // Dynamic structure
1077: typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1078: // Iterator interface
1079: typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1080: typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1081: // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1082: typedef std::set<point_type> pointSet;
1083: typedef ALE::array<point_type> pointArray;
1084: typedef pointSet coneSet;
1085: typedef pointSet supportSet;
1086: typedef pointArray coneArray;
1087: typedef pointArray supportArray;
1088: protected:
1089: // Arrow Containers
1090: typedef index_type *offsets_type;
1091: typedef point_type *cones_type;
1092: typedef point_type *supports_type;
1093: // Decorators
1094: typedef int *orientations_type;
1095: protected:
1096: // Data
1097: bool indexAllocated;
1098: offsets_type coneOffsets;
1099: offsets_type supportOffsets;
1100: bool pointAllocated;
1101: index_type maxConeSize;
1102: index_type maxSupportSize;
1103: index_type baseSize;
1104: index_type capSize;
1105: cones_type cones;
1106: supports_type supports;
1107: bool orientCones;
1108: orientations_type coneOrientations;
1109: chart_type chart;
1110: int_allocator_type intAlloc;
1111: index_allocator_type indexAlloc;
1112: point_allocator_type pointAlloc;
1113: newpoints_type newCones;
1114: newpoints_type newSupports;
1115: // Sequences
1116: Obj<coneSequence> coneSeq;
1117: Obj<supportSequence> supportSeq;
1118: protected: // Memory Management
1119: void createIndices() {
1120: this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1121: this->coneOffsets -= this->chart.min();
1122: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1123: this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1124: this->supportOffsets -= this->chart.min();
1125: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1126: this->indexAllocated = true;
1127: };
1128: void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1129: if (*coneOffsets) {
1130: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1131: *coneOffsets += chart.min();
1132: indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1133: *coneOffsets = NULL;
1134: }
1135: if (*supportOffsets) {
1136: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1137: *supportOffsets += chart.min();
1138: indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1139: *supportOffsets = NULL;
1140: }
1141: };
1142: void destroyIndices() {
1143: this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1144: this->indexAllocated = false;
1145: this->maxConeSize = -1;
1146: this->maxSupportSize = -1;
1147: this->baseSize = -1;
1148: this->capSize = -1;
1149: };
1150: void createPoints() {
1151: this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1152: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1153: this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1154: for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1155: if (orientCones) {
1156: this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1157: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1158: }
1159: this->pointAllocated = true;
1160: };
1161: void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1162: if (*cones) {
1163: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1164: pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1165: *cones = NULL;
1166: }
1167: if (*supports) {
1168: for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1169: pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1170: *supports = NULL;
1171: }
1172: if (*coneOrientations) {
1173: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1174: intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1175: *coneOrientations = NULL;
1176: }
1177: };
1178: void destroyPoints() {
1179: this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1180: this->pointAllocated = false;
1181: };
1182: void prefixSum(const offsets_type array) {
1183: for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1184: array[p] = array[p] + array[p-1];
1185: }
1186: };
1187: void calculateBaseAndCapSize() {
1188: this->baseSize = 0;
1189: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1190: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1191: ++this->baseSize;
1192: }
1193: }
1194: this->capSize = 0;
1195: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1196: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1197: ++this->capSize;
1198: }
1199: }
1200: };
1201: public:
1202: 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) {
1203: this->coneSeq = new coneSequence(NULL, 0, 0);
1204: this->supportSeq = new supportSequence(NULL, 0, 0);
1205: };
1206: 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) {
1207: this->setChart(chart_type(min, max));
1208: this->coneSeq = new coneSequence(NULL, 0, 0);
1209: this->supportSeq = new supportSequence(NULL, 0, 0);
1210: };
1211: ~IFSieve() {
1212: this->destroyPoints();
1213: this->destroyIndices();
1214: };
1215: public: // Accessors
1216: const chart_type& getChart() const {return this->chart;};
1217: void setChart(const chart_type& chart) {
1218: this->destroyPoints();
1219: this->destroyIndices();
1220: this->chart = chart;
1221: this->createIndices();
1222: };
1223: index_type getMaxConeSize() const {
1224: return this->maxConeSize;
1225: };
1226: index_type getMaxSupportSize() const {
1227: return this->maxSupportSize;
1228: };
1229: bool baseContains(const point_type& p) const {
1230: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1231: this->chart.checkPoint(p);
1233: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1234: return true;
1235: }
1236: return false;
1237: };
1238: bool orientedCones() const {return this->orientCones;};
1239: // Raw array access
1240: offsets_type getConeOffsets() {return this->coneOffsets;};
1241: offsets_type getSupportOffsets() {return this->supportOffsets;};
1242: cones_type getCones() {return this->cones;};
1243: supports_type getSupports() {return this->supports;};
1244: orientations_type getConeOrientations() {return this->coneOrientations;};
1245: public: // Construction
1246: index_type getConeSize(const point_type& p) const {
1247: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1248: this->chart.checkPoint(p);
1249: return this->coneOffsets[p+1]-this->coneOffsets[p];
1250: };
1251: void setConeSize(const point_type& p, const index_type c) {
1252: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1253: this->chart.checkPoint(p);
1254: this->coneOffsets[p+1] = c;
1255: this->maxConeSize = std::max(this->maxConeSize, c);
1256: };
1257: void addConeSize(const point_type& p, const index_type c) {
1258: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1259: this->chart.checkPoint(p);
1260: this->coneOffsets[p+1] += c;
1261: this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1262: };
1263: index_type getSupportSize(const point_type& p) const {
1264: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1265: this->chart.checkPoint(p);
1266: return this->supportOffsets[p+1]-this->supportOffsets[p];
1267: };
1268: void setSupportSize(const point_type& p, const index_type s) {
1269: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1270: this->chart.checkPoint(p);
1271: this->supportOffsets[p+1] = s;
1272: this->maxSupportSize = std::max(this->maxSupportSize, s);
1273: };
1274: void addSupportSize(const point_type& p, const index_type s) {
1275: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1276: this->chart.checkPoint(p);
1277: this->supportOffsets[p+1] += s;
1278: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1279: };
1280: void allocate() {
1281: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1282: this->prefixSum(this->coneOffsets);
1283: this->prefixSum(this->supportOffsets);
1284: this->createPoints();
1285: this->calculateBaseAndCapSize();
1286: }
1287: // Purely for backwards compatibility
1288: template<typename Color>
1289: void addArrow(const point_type& p, const point_type& q, const Color c, const bool forceSupport = false) {
1290: this->addArrow(p, q, forceSupport);
1291: }
1292: void addArrow(const point_type& p, const point_type& q, const bool forceSupport = false) {
1293: if (!this->chart.hasPoint(q)) {
1294: if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1295: const index_type start = this->coneOffsets[q];
1296: const index_type end = this->coneOffsets[q+1];
1298: for(int c = start; c < end; ++c) {
1299: this->newCones[q].push_back(this->cones[c]);
1300: }
1301: }
1302: this->newCones[q].push_back(p);
1303: }
1304: if (!this->chart.hasPoint(p) || forceSupport) {
1305: if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1306: const index_type start = this->supportOffsets[p];
1307: const index_type end = this->supportOffsets[p+1];
1309: for(int s = start; s < end; ++s) {
1310: this->newSupports[p].push_back(this->supports[s]);
1311: }
1312: }
1313: this->newSupports[p].push_back(q);
1314: }
1315: };
1316: void reallocate() {
1317: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1318: if (!this->newCones.size() && !this->newSupports.size()) return;
1319: const chart_type oldChart = this->chart;
1320: offsets_type oldConeOffsets = this->coneOffsets;
1321: offsets_type oldSupportOffsets = this->supportOffsets;
1322: cones_type oldCones = this->cones;
1323: supports_type oldSupports = this->supports;
1324: orientations_type oldConeOrientations = this->coneOrientations;
1325: point_type min = this->chart.min();
1326: point_type max = this->chart.max()-1;
1328: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1329: min = std::min(min, c_iter->first);
1330: max = std::max(max, c_iter->first);
1331: }
1332: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1333: min = std::min(min, s_iter->first);
1334: max = std::max(max, s_iter->first);
1335: }
1336: this->chart = chart_type(min, max+1);
1337: this->createIndices();
1338: // Copy sizes (converted from offsets)
1339: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1340: this->coneOffsets[p+1] = oldConeOffsets[p+1]-oldConeOffsets[p];
1341: this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1342: }
1343: // Inject new sizes
1344: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1345: this->coneOffsets[c_iter->first+1] = c_iter->second.size();
1346: this->maxConeSize = std::max(this->maxConeSize, (int) c_iter->second.size());
1347: }
1348: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1349: this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1350: this->maxSupportSize = std::max(this->maxSupportSize, (int) s_iter->second.size());
1351: }
1352: this->prefixSum(this->coneOffsets);
1353: this->prefixSum(this->supportOffsets);
1354: this->createPoints();
1355: this->calculateBaseAndCapSize();
1356: // Copy cones and supports
1357: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1358: const index_type cStart = this->coneOffsets[p];
1359: const index_type cOStart = oldConeOffsets[p];
1360: const index_type cOEnd = oldConeOffsets[p+1];
1361: const index_type sStart = this->supportOffsets[p];
1362: const index_type sOStart = oldSupportOffsets[p];
1363: const index_type sOEnd = oldSupportOffsets[p+1];
1365: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1366: this->cones[c] = oldCones[cO];
1367: }
1368: for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1369: this->supports[s] = oldSupports[sO];
1370: }
1371: if (this->orientCones) {
1372: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1373: this->coneOrientations[c] = oldConeOrientations[cO];
1374: }
1375: }
1376: }
1377: // Inject new cones and supports
1378: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1379: index_type start = this->coneOffsets[c_iter->first];
1381: for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1382: this->cones[start++] = *p_iter;
1383: }
1384: if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1385: }
1386: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1387: index_type start = this->supportOffsets[s_iter->first];
1389: for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1390: this->supports[start++] = *p_iter;
1391: }
1392: if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1393: }
1394: this->newCones.clear();
1395: this->newSupports.clear();
1396: this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1397: this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1398: };
1399: // Recalculate the support offsets and fill the supports
1400: // This is used when an IFSieve is being used as a label
1401: void recalculateLabel() {
1402: ISieveVisitor::PointRetriever<IFSieve> v(1);
1404: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1405: this->supportOffsets[p+1] = 0;
1406: }
1407: this->maxSupportSize = 0;
1408: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1409: this->cone(p, v);
1410: if (v.getSize()) {
1411: this->supportOffsets[v.getPoints()[0]+1]++;
1412: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1413: }
1414: v.clear();
1415: }
1416: this->prefixSum(this->supportOffsets);
1417: this->calculateBaseAndCapSize();
1418: this->symmetrize();
1419: };
1420: void setCone(const point_type cone[], const point_type& p) {
1421: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1422: this->chart.checkPoint(p);
1423: const index_type start = this->coneOffsets[p];
1424: const index_type end = this->coneOffsets[p+1];
1426: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1427: this->cones[c] = cone[i];
1428: }
1429: };
1430: void setCone(const point_type cone, const point_type& p) {
1431: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1432: this->chart.checkPoint(p);
1433: const index_type start = this->coneOffsets[p];
1434: const index_type end = this->coneOffsets[p+1];
1436: if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
1437: this->cones[start] = cone;
1438: };
1439: #if 0
1440: template<typename PointSequence>
1441: void setCone(const PointSequence& cone, const point_type& p) {
1442: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1443: this->chart.checkPoint(p);
1444: const index_type start = this->coneOffsets[p];
1445: const index_type end = this->coneOffsets[p+1];
1446: if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
1447: typename PointSequence::iterator c_iter = cone.begin();
1449: for(index_type c = start; c < end; ++c, ++c_iter) {
1450: this->cones[c] = *c_iter;
1451: }
1452: };
1453: #endif
1454: void setSupport(const point_type& p, const point_type support[]) {
1455: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1456: this->chart.checkPoint(p);
1457: const index_type start = this->supportOffsets[p];
1458: const index_type end = this->supportOffsets[p+1];
1460: for(index_type s = start, i = 0; s < end; ++s, ++i) {
1461: this->supports[s] = support[i];
1462: }
1463: };
1464: #if 0
1465: template<typename PointSequence>
1466: void setSupport(const point_type& p, const PointSequence& support) {
1467: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1468: this->chart.checkPoint(p);
1469: const index_type start = this->supportOffsets[p];
1470: const index_type end = this->supportOffsets[p+1];
1471: if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
1472: typename PointSequence::iterator s_iter = support.begin();
1474: for(index_type s = start; s < end; ++s, ++s_iter) {
1475: this->supports[s] = *s_iter;
1476: }
1477: };
1478: #endif
1479: void setConeOrientation(const int coneOrientation[], const point_type& p) {
1480: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1481: this->chart.checkPoint(p);
1482: const index_type start = this->coneOffsets[p];
1483: const index_type end = this->coneOffsets[p+1];
1485: for(index_type c = start, i = 0; c < end; ++c, ++i) {
1486: this->coneOrientations[c] = coneOrientation[i];
1487: }
1488: };
1489: void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
1490: for(point_type p = 0; p < numCells; ++p) {
1491: const index_type start = p*numCorners;
1492: const index_type end = (p+1)*numCorners;
1494: for(index_type c = start; c < end; ++c) {
1495: const point_type q = cones[c]+offset;
1497: this->supportOffsets[q+1]++;
1498: }
1499: }
1500: for(point_type p = numCells; p < this->chart.max(); ++p) {
1501: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1502: }
1503: };
1504: void symmetrize() {
1505: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1506: offsets -= this->chart.min();
1507: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1508: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1510: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1511: const index_type start = this->coneOffsets[p];
1512: const index_type end = this->coneOffsets[p+1];
1514: for(index_type c = start; c < end; ++c) {
1515: const point_type q = this->cones[c];
1517: this->chart.checkPoint(q);
1518: this->supports[this->supportOffsets[q]+offsets[q]] = p;
1519: ++offsets[q];
1520: }
1521: }
1522: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1523: indexAlloc.deallocate(offsets, this->chart.size()+1);
1524: }
1525: index_type getBaseSize() const {
1526: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1527: return this->baseSize;
1528: }
1529: index_type getCapSize() const {
1530: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1531: return this->capSize;
1532: }
1533: template<typename _Section>
1534: void relabel(_Section& labeling) {
1537: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
1538: offsets -= this->chart.min();
1539: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
1540: // Recalculate coneOffsets
1541: for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
1542: const point_type newP = labeling.restrictPoint(p)[0];
1544: offsets[newP+1] = this->getConeSize(p);
1545: }
1546: this->prefixSum(offsets);
1547: PetscMemcpy(this->coneOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
1548: // Recalculate supportOffsets
1549: for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
1550: const point_type newP = labeling.restrictPoint(p)[0];
1552: offsets[newP+1] = this->getSupportSize(p);
1553: }
1554: this->prefixSum(offsets);
1555: PetscMemcpy(this->supportOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
1556: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
1557: indexAlloc.deallocate(offsets, this->chart.size()+1);
1558: index_type size = std::max(this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()],
1559: this->supportOffsets[this->chart.max()] - this->supportOffsets[this->chart.min()]);
1560: index_type *orientations = offsets = indexAlloc.allocate(size);
1561: for(index_type i = 0; i < size; ++i) {indexAlloc.construct(orientations+i, index_type(0));}
1562: // Recalculate coneOrientations
1563: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1564: const point_type newP = labeling.restrictPoint(p)[0];
1565: const index_type start = this->coneOffsets[newP];
1566: const index_type end = this->coneOffsets[newP+1];
1568: for(index_type c = start; c < end; ++c, ++offset) {
1569: orientations[c] = this->coneOrientations[offset];
1570: }
1571: }
1572: PetscMemcpy(this->coneOrientations, orientations, (this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()])*sizeof(index_type));CHKERRXX(ierr);
1573: for(index_type i = 0; i < size; ++i) {indexAlloc.destroy(orientations+i);}
1574: indexAlloc.deallocate(orientations, size);
1575: // Recalculate cones
1576: point_type *array = pointAlloc.allocate(size);
1578: for(index_type i = 0; i < size; ++i) {pointAlloc.construct(array+i, point_type(0));}
1579: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1580: const point_type newP = labeling.restrictPoint(p)[0];
1581: const index_type start = this->coneOffsets[newP];
1582: const index_type end = this->coneOffsets[newP+1];
1584: for(index_type c = start; c < end; ++c, ++offset) {
1585: const point_type newQ = labeling.restrictPoint(this->cones[offset])[0];
1587: array[c] = newQ;
1588: }
1589: }
1590: PetscMemcpy(this->cones, array, size*sizeof(point_type));CHKERRXX(ierr);
1591: // Recalculate supports
1592: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
1593: const point_type newP = labeling.restrictPoint(p)[0];
1594: const index_type start = this->supportOffsets[newP];
1595: const index_type end = this->supportOffsets[newP+1];
1597: for(index_type c = start; c < end; ++c, ++offset) {
1598: const point_type newQ = labeling.restrictPoint(this->supports[offset])[0];
1600: array[c] = newQ;
1601: }
1602: }
1603: PetscMemcpy(this->supports, array, size*sizeof(point_type));CHKERRXX(ierr);
1604: for(index_type i = 0; i < size; ++i) {pointAlloc.destroy(array+i);}
1605: pointAlloc.deallocate(array, size);
1606: }
1607: public: // Traversals
1608: template<typename Visitor>
1609: void roots(const Visitor& v) const {
1610: this->roots(const_cast<Visitor&>(v));
1611: }
1612: template<typename Visitor>
1613: void roots(Visitor& v) const {
1614: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1616: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1617: if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
1618: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1619: v.visitPoint(p);
1620: }
1621: }
1622: }
1623: }
1624: template<typename Visitor>
1625: void leaves(const Visitor& v) const {
1626: this->leaves(const_cast<Visitor&>(v));
1627: }
1628: template<typename Visitor>
1629: void leaves(Visitor& v) const {
1630: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1632: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1633: if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
1634: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1635: v.visitPoint(p);
1636: }
1637: }
1638: }
1639: }
1640: template<typename Visitor>
1641: void base(const Visitor& v) const {
1642: this->base(const_cast<Visitor&>(v));
1643: }
1644: template<typename Visitor>
1645: void base(Visitor& v) const {
1646: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1648: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1649: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1650: v.visitPoint(p);
1651: }
1652: }
1653: }
1654: template<typename Visitor>
1655: void cap(const Visitor& v) const {
1656: this->cap(const_cast<Visitor&>(v));
1657: }
1658: template<typename Visitor>
1659: void cap(Visitor& v) const {
1660: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1662: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1663: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1664: v.visitPoint(p);
1665: }
1666: }
1667: }
1668: template<typename PointSequence, typename Visitor>
1669: void cone(const PointSequence& points, Visitor& v) const {
1670: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1671: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1672: const point_type p = *p_iter;
1673: this->chart.checkPoint(p);
1674: const index_type start = this->coneOffsets[p];
1675: const index_type end = this->coneOffsets[p+1];
1677: for(index_type c = start; c < end; ++c) {
1678: v.visitArrow(arrow_type(this->cones[c], p));
1679: v.visitPoint(this->cones[c]);
1680: }
1681: }
1682: }
1683: template<typename Visitor>
1684: void cone(const point_type& p, Visitor& v) const {
1685: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1686: this->chart.checkPoint(p);
1687: const index_type start = this->coneOffsets[p];
1688: const index_type end = this->coneOffsets[p+1];
1690: for(index_type c = start; c < end; ++c) {
1691: v.visitArrow(arrow_type(this->cones[c], p));
1692: v.visitPoint(this->cones[c]);
1693: }
1694: }
1695: const Obj<coneSequence>& cone(const point_type& p) const {
1696: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1697: if (!this->chart.hasPoint(p)) {
1698: this->coneSeq->reset(this->cones, 0, 0);
1699: } else {
1700: this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
1701: }
1702: return this->coneSeq;
1703: }
1704: template<typename PointSequence, typename Visitor>
1705: void support(const PointSequence& points, Visitor& v) const {
1706: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1707: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1708: const point_type p = *p_iter;
1709: this->chart.checkPoint(p);
1710: const index_type start = this->supportOffsets[p];
1711: const index_type end = this->supportOffsets[p+1];
1713: for(index_type s = start; s < end; ++s) {
1714: v.visitArrow(arrow_type(p, this->supports[s]));
1715: v.visitPoint(this->supports[s]);
1716: }
1717: }
1718: }
1719: template<typename Visitor>
1720: void support(const point_type& p, Visitor& v) const {
1721: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1722: this->chart.checkPoint(p);
1723: const index_type start = this->supportOffsets[p];
1724: const index_type end = this->supportOffsets[p+1];
1726: for(index_type s = start; s < end; ++s) {
1727: v.visitArrow(arrow_type(p, this->supports[s]));
1728: v.visitPoint(this->supports[s]);
1729: }
1730: }
1731: const Obj<supportSequence>& support(const point_type& p) const {
1732: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1733: if (!this->chart.hasPoint(p)) {
1734: this->supportSeq->reset(this->supports, 0, 0);
1735: } else {
1736: this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
1737: }
1738: return this->supportSeq;
1739: }
1740: template<typename Visitor>
1741: void orientedCone(const point_type& p, Visitor& v) const {
1742: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1743: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1744: this->chart.checkPoint(p);
1745: const index_type start = this->coneOffsets[p];
1746: const index_type end = this->coneOffsets[p+1];
1748: for(index_type c = start; c < end; ++c) {
1749: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1750: v.visitPoint(this->cones[c], this->coneOrientations[c]);
1751: }
1752: }
1753: template<typename Visitor>
1754: void orientedReverseCone(const point_type& p, Visitor& v) const {
1755: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1756: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1757: this->chart.checkPoint(p);
1758: const index_type start = this->coneOffsets[p];
1759: const index_type end = this->coneOffsets[p+1];
1761: for(index_type c = end-1; c >= start; --c) {
1762: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
1763: v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
1764: }
1765: }
1766: template<typename Visitor>
1767: void orientedSupport(const point_type& p, Visitor& v) const {
1768: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1769: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1770: this->chart.checkPoint(p);
1771: const index_type start = this->supportOffsets[p];
1772: const index_type end = this->supportOffsets[p+1];
1774: for(index_type s = start; s < end; ++s) {
1775: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1776: //v.visitPoint(this->supports[s], this->supportOrientations[s]);
1777: v.visitArrow(arrow_type(this->supports[s], p), 0);
1778: v.visitPoint(this->supports[s], 0);
1779: }
1780: }
1781: template<typename Visitor>
1782: void orientedReverseSupport(const point_type& p, Visitor& v) const {
1783: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
1784: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1785: this->chart.checkPoint(p);
1786: const index_type start = this->supportOffsets[p];
1787: const index_type end = this->supportOffsets[p+1];
1789: for(index_type s = end-1; s >= start; --s) {
1790: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
1791: //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
1792: v.visitArrow(arrow_type(this->supports[s], p), 0);
1793: v.visitPoint(this->supports[s], 0);
1794: }
1795: }
1796: // Currently does only 1 level
1797: // Does not check for uniqueness
1798: template<typename Visitor>
1799: void meet(const point_type& p, const point_type& q, Visitor& v) const {
1800: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1801: this->chart.checkPoint(p);
1802: this->chart.checkPoint(q);
1803: const index_type startP = this->coneOffsets[p];
1804: const index_type endP = this->coneOffsets[p+1];
1805: const index_type startQ = this->coneOffsets[q];
1806: const index_type endQ = this->coneOffsets[q+1];
1808: for(index_type cP = startP; cP < endP; ++cP) {
1809: const point_type& c1 = this->cones[cP];
1811: for(index_type cQ = startQ; cQ < endQ; ++cQ) {
1812: if (c1 == this->cones[cQ]) v.visitPoint(c1);
1813: }
1814: if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
1815: }
1816: }
1817: // Currently does only 1 level
1818: template<typename Sequence, typename Visitor>
1819: void join(const Sequence& points, Visitor& v) const {
1820: typedef std::set<point_type> pointSet;
1821: pointSet intersect[2] = {pointSet(), pointSet()};
1822: pointSet tmp;
1823: int p = 0;
1824: int c = 0;
1826: for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
1827: this->chart.checkPoint(*p_iter);
1828: tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
1829: if (p == 0) {
1830: intersect[1-c].insert(tmp.begin(), tmp.end());
1831: p++;
1832: } else {
1833: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1834: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1835: intersect[c].clear();
1836: }
1837: c = 1 - c;
1838: tmp.clear();
1839: }
1840: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1841: v.visitPoint(*p_iter);
1842: }
1843: }
1844: // Helper function
1845: void insertNSupport(point_type p, pointSet& set, const int depth) {
1846: const index_type start = this->supportOffsets[p];
1847: const index_type end = this->supportOffsets[p+1];
1849: if (depth == 1) {
1850: set.insert(&this->supports[start], &this->supports[end]);
1851: } else {
1852: for(index_type s = start; s < end; ++s) {
1853: this->insertNSupport(this->supports[s], set, depth-1);
1854: }
1855: }
1856: }
1857: // Gives only the join of depth n
1858: template<typename SequenceIterator, typename Visitor>
1859: void nJoin(const SequenceIterator& pointsBegin, const SequenceIterator& pointsEnd, const int depth, Visitor& v) {
1860: typedef std::set<point_type> pointSet;
1861: pointSet intersect[2] = {pointSet(), pointSet()};
1862: pointSet tmp;
1863: int p = 0;
1864: int c = 0;
1866: for(SequenceIterator p_iter = pointsBegin; p_iter != pointsEnd; ++p_iter) {
1867: this->chart.checkPoint(*p_iter);
1868: // Put points in the nSupport into tmp (duplicates are fine since it is a set)
1869: this->insertNSupport(*p_iter, tmp, depth);
1870: if (p == 0) {
1871: intersect[1-c].insert(tmp.begin(), tmp.end());
1872: p++;
1873: } else {
1874: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
1875: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
1876: intersect[c].clear();
1877: }
1878: c = 1 - c;
1879: tmp.clear();
1880: }
1881: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
1882: v.visitPoint(*p_iter);
1883: }
1884: }
1885: public: // Viewing
1886: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
1887: ostringstream txt;
1888: int rank;
1890: if (comm == MPI_COMM_NULL) {
1891: comm = this->comm();
1892: rank = this->commRank();
1893: } else {
1894: MPI_Comm_rank(comm, &rank);
1895: }
1896: if (name == "") {
1897: if(rank == 0) {
1898: txt << "viewing an IFSieve" << std::endl;
1899: }
1900: } else {
1901: if(rank == 0) {
1902: txt << "viewing IFSieve '" << name << "'" << std::endl;
1903: }
1904: }
1905: PetscSynchronizedPrintf(comm, "Max sizes cone: %d support: %d\n", this->getMaxConeSize(), this->getMaxSupportSize());
1906: if(rank == 0) {
1907: txt << "cap --> base:" << std::endl;
1908: }
1909: ISieveVisitor::PrintVisitor pV(txt, rank);
1910: this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
1911: PetscSynchronizedPrintf(comm, txt.str().c_str());
1912: PetscSynchronizedFlush(comm);
1913: ostringstream txt2;
1915: if(rank == 0) {
1916: txt2 << "base <-- cap:" << std::endl;
1917: }
1918: ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
1919: this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
1920: PetscSynchronizedPrintf(comm, txt2.str().c_str());
1921: PetscSynchronizedFlush(comm);
1922: if (orientCones) {
1923: ostringstream txt3;
1925: if(rank == 0) {
1926: txt3 << "Orientation:" << std::endl;
1927: }
1928: ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
1929: this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
1930: PetscSynchronizedPrintf(comm, txt3.str().c_str());
1931: PetscSynchronizedFlush(comm);
1932: }
1933: };
1934: };
1936: class ISieveConverter {
1937: public:
1938: template<typename Sieve, typename ISieve, typename Renumbering>
1939: static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
1940: // First construct a renumbering of the sieve points
1941: const Obj<typename Sieve::baseSequence>& base = sieve.base();
1942: const Obj<typename Sieve::capSequence>& cap = sieve.cap();
1943: typename ISieve::point_type min = 0;
1944: typename ISieve::point_type max = 0;
1946: if (renumber) {
1947: /* Roots/Leaves from Sieve do not seem to work */
1949: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1950: if (sieve.support(*b_iter)->size() == 0) {
1951: renumbering[*b_iter] = max++;
1952: }
1953: }
1954: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1955: if (sieve.cone(*c_iter)->size() == 0) {
1956: renumbering[*c_iter] = max++;
1957: }
1958: }
1959: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1960: if (renumbering.find(*b_iter) == renumbering.end()) {
1961: renumbering[*b_iter] = max++;
1962: }
1963: }
1964: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1965: if (renumbering.find(*c_iter) == renumbering.end()) {
1966: renumbering[*c_iter] = max++;
1967: }
1968: }
1969: } else {
1970: if (base->size()) {
1971: min = *base->begin();
1972: max = *base->begin();
1973: } else if (cap->size()) {
1974: min = *cap->begin();
1975: max = *cap->begin();
1976: }
1977: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1978: min = std::min(min, *b_iter);
1979: max = std::max(max, *b_iter);
1980: }
1981: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
1982: min = std::min(min, *c_iter);
1983: max = std::max(max, *c_iter);
1984: }
1985: if (base->size() || cap->size()) {
1986: ++max;
1987: }
1988: for(typename ISieve::point_type p = min; p < max; ++p) {
1989: renumbering[p] = p;
1990: }
1991: }
1992: // Create the ISieve
1993: isieve.setChart(typename ISieve::chart_type(min, max));
1994: // Set cone and support sizes
1995: size_t maxSize = 0;
1997: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1998: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2000: isieve.setConeSize(renumbering[*b_iter], cone->size());
2001: maxSize = std::max(maxSize, cone->size());
2002: }
2003: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2004: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2006: isieve.setSupportSize(renumbering[*c_iter], support->size());
2007: maxSize = std::max(maxSize, support->size());
2008: }
2009: isieve.allocate();
2010: // Fill up cones and supports
2011: typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];
2013: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2014: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2015: int i = 0;
2017: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2018: points[i] = renumbering[*c_iter];
2019: }
2020: isieve.setCone(points, renumbering[*b_iter]);
2021: }
2022: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2023: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2024: int i = 0;
2026: for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
2027: points[i] = renumbering[*s_iter];
2028: }
2029: isieve.setSupport(renumbering[*c_iter], points);
2030: }
2031: delete [] points;
2032: }
2033: template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
2034: static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
2035: if (isieve.getMaxConeSize() < 0) return;
2036: const Obj<typename Sieve::baseSequence>& base = sieve.base();
2037: int *orientations = new int[isieve.getMaxConeSize()];
2039: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2040: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2041: int i = 0;
2043: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2044: typename ArrowSection::point_type arrow(*c_iter, *b_iter);
2046: orientations[i] = orientation->restrictPoint(arrow)[0];
2047: }
2048: isieve.setConeOrientation(orientations, renumbering[*b_iter]);
2049: }
2050: delete [] orientations;
2051: }
2052: template<typename Section, typename ISection, typename Renumbering>
2053: static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
2054: const typename Section::chart_type& chart = coordinates.getChart();
2055: typename ISection::point_type min = *chart.begin();
2056: typename ISection::point_type max = *chart.begin();
2058: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2059: min = std::min(min, renumbering[*p_iter]);
2060: max = std::max(max, renumbering[*p_iter]);
2061: }
2062: icoordinates.setChart(typename ISection::chart_type(min, max+1));
2063: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2064: icoordinates.setFiberDimension(renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));
2065: }
2066: icoordinates.allocatePoint();
2067: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2068: icoordinates.updatePoint(renumbering[*p_iter], coordinates.restrictPoint(*p_iter));
2069: }
2070: }
2071: template<typename Label, typename Renumbering>
2072: static void convertLabel(const Obj<Label>& newLabel, const Obj<Label>& oldLabel, Renumbering& renumbering) {
2073: for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
2074: if (oldLabel->getConeSize(p->first)) {
2075: newLabel->setCone(*oldLabel->cone(p->first)->begin(), p->second);
2076: }
2077: }
2078: }
2079: template<typename Mesh, typename IMesh, typename Renumbering>
2080: static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
2081: convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
2082: imesh.stratify();
2083: convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
2084: convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
2085: if (mesh.hasRealSection("normals")) {
2086: convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
2087: }
2088: const typename Mesh::labels_type& labels = mesh.getLabels();
2089: std::string heightName("height");
2090: std::string depthName("depth");
2092: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
2093: #ifdef IMESH_NEW_LABELS
2094: if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
2095: convertLabel(imesh, l_iter->first, l_iter->second);
2096: }
2097: #else
2098: if (renumber) {
2099: convertLabel(imesh.createLabel(l_iter->first), l_iter->second, renumbering);
2100: } else {
2101: imesh.setLabel(l_iter->first, l_iter->second);
2102: }
2103: #endif
2104: }
2105: }
2106: };
2108: class ISieveSerializer {
2109: public:
2110: template<typename ISieve>
2111: static void writeSieve(const std::string& filename, ISieve& sieve) {
2112: std::ofstream fs;
2114: if (sieve.commRank() == 0) {
2115: fs.open(filename.c_str());
2116: }
2117: writeSieve(fs, sieve);
2118: if (sieve.commRank() == 0) {
2119: fs.close();
2120: }
2121: }
2122: template<typename ISieve>
2123: static void writeSieve(std::ofstream& fs, ISieve& sieve) {
2124: typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
2125: const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
2126: typename ISieve::point_type min = chart->min();
2127: typename ISieve::point_type max = chart->max();
2128: PetscInt *mins = new PetscInt[sieve.commSize()];
2129: PetscInt *maxs = new PetscInt[sieve.commSize()];
2131: // Write sizes
2132: if (sieve.commRank() == 0) {
2133: // Write local
2134: fs << min <<" "<< max << std::endl;
2135: for(typename ISieve::point_type p = min; p < max; ++p) {
2136: fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
2137: }
2138: // Receive and write remote
2139: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2140: PetscInt minp, maxp;
2141: PetscInt *sizes;
2142: MPI_Status status;
2145: MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2146: MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2147: PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2148: MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2149: fs << minp <<" "<< maxp << std::endl;
2150: for(PetscInt s = 0; s < maxp - minp; ++s) {
2151: fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
2152: }
2153: PetscFree(sizes);CHKERRXX(ierr);
2154: mins[pr] = minp;
2155: maxs[pr] = maxp;
2156: }
2157: } else {
2158: // Send remote
2159: //PetscInt min = chart->min();
2160: //PetscInt max = chart->max();
2161: PetscInt s = 0;
2162: PetscInt *sizes;
2165: MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2166: MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2167: PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
2168: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2169: sizes[s*2+0] = sieve.getConeSize(p);
2170: sizes[s*2+1] = sieve.getSupportSize(p);
2171: }
2172: MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2173: PetscFree(sizes);CHKERRXX(ierr);
2174: }
2175: // Write covering
2176: if (sieve.commRank() == 0) {
2177: // Write local
2178: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2180: for(typename ISieve::point_type p = min; p < max; ++p) {
2181: sieve.cone(p, pV);
2182: const typename Visitor::point_type *cone = pV.getPoints();
2183: const int cSize = pV.getSize();
2185: if (cSize > 0) {
2186: for(int c = 0; c < cSize; ++c) {
2187: if (c) {fs << " ";}
2188: fs << cone[c];
2189: }
2190: fs << std::endl;
2191: }
2192: pV.clear();
2194: sieve.orientedCone(p, pV);
2195: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2196: const int oSize = pV.getOrientedSize();
2198: if (oSize > 0) {
2199: for(int c = 0; c < oSize; ++c) {
2200: if (c) {fs << " ";}
2201: fs << oCone[c].second;
2202: }
2203: fs << std::endl;
2204: }
2205: pV.clear();
2207: sieve.support(p, pV);
2208: const typename Visitor::point_type *support = pV.getPoints();
2209: const int sSize = pV.getSize();
2211: if (sSize > 0) {
2212: for(int s = 0; s < sSize; ++s) {
2213: if (s) {fs << " ";}
2214: fs << support[s];
2215: }
2216: fs << std::endl;
2217: }
2218: pV.clear();
2219: }
2220: // Receive and write remote
2221: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2222: PetscInt size;
2223: PetscInt *data;
2224: PetscInt off = 0;
2225: MPI_Status status;
2228: MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2229: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2230: MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2231: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2232: PetscInt cSize = data[off++];
2234: fs << cSize << std::endl;
2235: if (cSize > 0) {
2236: for(int c = 0; c < cSize; ++c) {
2237: if (c) {fs << " ";}
2238: fs << data[off++];
2239: }
2240: fs << std::endl;
2241: }
2242: PetscInt oSize = data[off++];
2244: if (oSize > 0) {
2245: for(int c = 0; c < oSize; ++c) {
2246: if (c) {fs << " ";}
2247: fs << data[off++];
2248: }
2249: fs << std::endl;
2250: }
2251: PetscInt sSize = data[off++];
2253: fs << sSize << std::endl;
2254: if (sSize > 0) {
2255: for(int s = 0; s < sSize; ++s) {
2256: if (s) {fs << " ";}
2257: fs << data[off++];
2258: }
2259: fs << std::endl;
2260: }
2261: }
2262: assert(off == size);
2263: PetscFree(data);CHKERRXX(ierr);
2264: }
2265: } else {
2266: // Send remote
2267: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
2268: PetscInt totalConeSize = 0;
2269: PetscInt totalSupportSize = 0;
2271: for(typename ISieve::point_type p = min; p < max; ++p) {
2272: totalConeSize += sieve.getConeSize(p);
2273: totalSupportSize += sieve.getSupportSize(p);
2274: }
2275: PetscInt size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
2276: PetscInt off = 0;
2277: PetscInt *data;
2280: MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2281: // There is no nice way to make a generic MPI type here. Really sucky
2282: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2283: for(typename ISieve::point_type p = min; p < max; ++p) {
2284: sieve.cone(p, pV);
2285: const typename Visitor::point_type *cone = pV.getPoints();
2286: const int cSize = pV.getSize();
2288: data[off++] = cSize;
2289: for(int c = 0; c < cSize; ++c) {
2290: data[off++] = cone[c];
2291: }
2292: pV.clear();
2294: sieve.orientedCone(p, pV);
2295: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
2296: const int oSize = pV.getOrientedSize();
2298: data[off++] = oSize;
2299: for(int c = 0; c < oSize; ++c) {
2300: data[off++] = oCone[c].second;
2301: }
2302: pV.clear();
2304: sieve.support(p, pV);
2305: const typename Visitor::point_type *support = pV.getPoints();
2306: const int sSize = pV.getSize();
2308: data[off++] = sSize;
2309: for(int s = 0; s < sSize; ++s) {
2310: data[off++] = support[s];
2311: }
2312: pV.clear();
2313: }
2314: MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
2315: PetscFree(data);CHKERRXX(ierr);
2316: }
2317: delete [] mins;
2318: delete [] maxs;
2319: // Output renumbering
2320: }
2321: template<typename ISieve>
2322: static void loadSieve(const std::string& filename, ISieve& sieve) {
2323: std::ifstream fs;
2325: if (sieve.commRank() == 0) {
2326: fs.open(filename.c_str());
2327: }
2328: loadSieve(fs, sieve);
2329: if (sieve.commRank() == 0) {
2330: fs.close();
2331: }
2332: }
2333: template<typename ISieve>
2334: static void loadSieve(std::ifstream& fs, ISieve& sieve) {
2335: typename ISieve::point_type min, max;
2336: PetscInt *mins = new PetscInt[sieve.commSize()];
2337: PetscInt *maxs = new PetscInt[sieve.commSize()];
2338: PetscInt *totalConeSizes = new PetscInt[sieve.commSize()];
2339: PetscInt *totalSupportSizes = new PetscInt[sieve.commSize()];
2341: // Load sizes
2342: if (sieve.commRank() == 0) {
2343: // Load local
2344: fs >> min;
2345: fs >> max;
2346: sieve.setChart(typename ISieve::chart_type(min, max));
2347: for(typename ISieve::point_type p = min; p < max; ++p) {
2348: typename ISieve::index_type coneSize, supportSize;
2350: fs >> coneSize;
2351: fs >> supportSize;
2352: sieve.setConeSize(p, coneSize);
2353: sieve.setSupportSize(p, supportSize);
2354: }
2355: // Load and send remote
2356: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2357: PetscInt minp, maxp;
2358: PetscInt *sizes;
2361: fs >> minp;
2362: fs >> maxp;
2363: MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2364: MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2365: PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2366: totalConeSizes[pr] = 0;
2367: totalSupportSizes[pr] = 0;
2368: for(PetscInt s = 0; s < maxp - minp; ++s) {
2369: fs >> sizes[s*2+0];
2370: fs >> sizes[s*2+1];
2371: totalConeSizes[pr] += sizes[s*2+0];
2372: totalSupportSizes[pr] += sizes[s*2+1];
2373: }
2374: MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2375: PetscFree(sizes);CHKERRXX(ierr);
2376: mins[pr] = minp;
2377: maxs[pr] = maxp;
2378: }
2379: } else {
2380: // Load remote
2381: PetscInt s = 0;
2382: PetscInt *sizes;
2383: MPI_Status status;
2386: MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2387: MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2388: sieve.setChart(typename ISieve::chart_type(min, max));
2389: PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2390: MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2391: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
2392: sieve.setConeSize(p, sizes[s*2+0]);
2393: sieve.setSupportSize(p, sizes[s*2+1]);
2394: }
2395: PetscFree(sizes);CHKERRXX(ierr);
2396: }
2397: sieve.allocate();
2398: // Load covering
2399: if (sieve.commRank() == 0) {
2400: // Load local
2401: typename ISieve::index_type maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
2402: typename ISieve::point_type *points = new typename ISieve::point_type[maxSize];
2404: for(typename ISieve::point_type p = min; p < max; ++p) {
2405: typename ISieve::index_type coneSize = sieve.getConeSize(p);
2406: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2408: if (coneSize > 0) {
2409: for(int c = 0; c < coneSize; ++c) {
2410: fs >> points[c];
2411: }
2412: sieve.setCone(points, p);
2413: if (sieve.orientedCones()) {
2414: for(int c = 0; c < coneSize; ++c) {
2415: fs >> points[c];
2416: }
2417: sieve.setConeOrientation(points, p);
2418: }
2419: }
2420: if (supportSize > 0) {
2421: for(int s = 0; s < supportSize; ++s) {
2422: fs >> points[s];
2423: }
2424: sieve.setSupport(p, points);
2425: }
2426: }
2427: delete [] points;
2428: // Load and send remote
2429: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2430: PetscInt size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
2431: PetscInt off = 0;
2432: PetscInt *data;
2435: MPI_Send(&size, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2436: // There is no nice way to make a generic MPI type here. Really sucky
2437: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
2438: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
2439: PetscInt coneSize, supportSize;
2441: fs >> coneSize;
2442: data[off++] = coneSize;
2443: if (coneSize > 0) {
2444: for(int c = 0; c < coneSize; ++c) {
2445: fs >> data[off++];
2446: }
2447: for(int c = 0; c < coneSize; ++c) {
2448: fs >> data[off++];
2449: }
2450: }
2451: fs >> supportSize;
2452: data[off++] = supportSize;
2453: if (supportSize > 0) {
2454: for(int s = 0; s < supportSize; ++s) {
2455: fs >> data[off++];
2456: }
2457: }
2458: }
2459: assert(off == size);
2460: MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
2461: PetscFree(data);CHKERRXX(ierr);
2462: }
2463: delete [] mins;
2464: delete [] maxs;
2465: } else {
2466: // Load remote
2467: PetscInt size;
2468: PetscInt *data;
2469: PetscInt off = 0;
2470: MPI_Status status;
2473: MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2474: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
2475: MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
2476: for(typename ISieve::point_type p = min; p < max; ++p) {
2477: typename ISieve::index_type coneSize = sieve.getConeSize(p);
2478: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
2479: PetscInt cs = data[off++];
2481: assert(cs == coneSize);
2482: if (coneSize > 0) {
2483: sieve.setCone(&data[off], p);
2484: off += coneSize;
2485: if (sieve.orientedCones()) {
2486: sieve.setConeOrientation(&data[off], p);
2487: off += coneSize;
2488: }
2489: }
2490: PetscInt ss = data[off++];
2492: assert(ss == supportSize);
2493: if (supportSize > 0) {
2494: sieve.setSupport(p, &data[off]);
2495: off += supportSize;
2496: }
2497: }
2498: assert(off == size);
2499: }
2500: // Load renumbering
2501: }
2502: };
2503: }
2505: #endif