Actual source code: ISieve.hh
petsc-3.4.5 2014-06-29
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 <petscdmplex.h>
10: #include <fstream>
12: //#define IMESH_NEW_LABELS
14: namespace ALE {
15: template<typename Point>
16: class OrientedPoint : public std::pair<Point, int> {
17: public:
18: OrientedPoint(const int o) : std::pair<Point, int>(o, o) {};
19: ~OrientedPoint() {};
20: friend std::ostream& operator<<(std::ostream& stream, const OrientedPoint& point) {
21: stream << "(" << point.first << ", " << point.second << ")";
22: return stream;
23: };
24: };
26: template<typename Point_, typename Alloc_ = malloc_allocator<Point_> >
27: class Interval {
28: public:
29: typedef Point_ point_type;
30: typedef Alloc_ alloc_type;
31: public:
32: class const_iterator {
33: protected:
34: point_type _p;
35: public:
36: const_iterator(const point_type p): _p(p) {};
37: ~const_iterator() {};
38: public:
39: const_iterator& operator=(const const_iterator& iter) {this->_p = iter._p; return *this;};
40: bool operator==(const const_iterator& iter) const {return this->_p == iter._p;};
41: bool operator!=(const const_iterator& iter) const {return this->_p != iter._p;};
42: const_iterator& operator++() {++this->_p; return *this;}
43: const_iterator& operator++(int) {
44: const_iterator tmp(*this);
45: ++(*this);
46: return tmp;
47: };
48: const_iterator& operator--() {--this->_p; return *this;}
49: const_iterator& operator--(int) {
50: const_iterator tmp(*this);
51: --(*this);
52: return tmp;
53: };
54: point_type operator*() const {return this->_p;};
55: };
56: protected:
57: point_type _min, _max;
58: public:
59: Interval(): _min(point_type()), _max(point_type()) {};
60: Interval(const point_type& min, const point_type& max): _min(min), _max(max) {};
61: Interval(const Interval& interval): _min(interval.min()), _max(interval.max()) {};
62: template<typename Iterator>
63: Interval(Iterator& iterator) {
64: this->_min = *std::min_element(iterator.begin(), iterator.end());
65: this->_max = (*std::max_element(iterator.begin(), iterator.end()))+1;
66: }
67: public:
68: Interval& operator=(const Interval& interval) {_min = interval.min(); _max = interval.max(); return *this;}
69: friend std::ostream& operator<<(std::ostream& stream, const Interval& interval) {
70: stream << "(" << interval.min() << ", " << interval.max() << ")";
71: return stream;
72: }
73: public:
74: const_iterator begin() const {return const_iterator(this->_min);};
75: const_iterator end() const {return const_iterator(this->_max);};
76: size_t size() const {return this->_max - this->_min;};
77: size_t count(const point_type& p) const {return ((p >= _min) && (p < _max)) ? 1 : 0;};
78: point_type min() const {return this->_min;};
79: point_type max() const {return this->_max;};
80: bool hasPoint(const point_type& point) const {
81: if (point < this->_min || point >= this->_max) return false;
82: return true;
83: };
84: void checkPoint(const point_type& point) const {
85: if (point < this->_min || point >= this->_max) {
86: ostringstream msg;
87: msg << "Invalid point " << point << " not in " << *this << std::endl;
88: throw ALE::Exception(msg.str().c_str());
89: }
90: };
91: };
93: template<typename Source_, typename Target_>
94: struct SimpleArrow {
95: typedef Source_ source_type;
96: typedef Target_ target_type;
97: const source_type source;
98: const target_type target;
99: SimpleArrow(const source_type& s, const target_type& t) : source(s), target(t) {};
100: template<typename OtherSource_, typename OtherTarget_>
101: struct rebind {
102: typedef SimpleArrow<OtherSource_, OtherTarget_> other;
103: };
104: struct flip {
105: typedef SimpleArrow<target_type, source_type> other;
106: other arrow(const SimpleArrow& a) {return type(a.target, a.source);};
107: };
108: friend bool operator<(const SimpleArrow& x, const SimpleArrow& y) {
109: return ((x.source < y.source) || ((x.source == y.source) && (x.target < y.target)));
110: };
111: friend std::ostream& operator<<(std::ostream& os, const SimpleArrow& a) {
112: os << a.source << " ----> " << a.target;
113: return os;
114: }
115: };
117: namespace ISieveVisitor {
118: template<typename Sieve>
119: class NullVisitor {
120: public:
121: inline void visitArrow(const typename Sieve::arrow_type&) {};
122: inline void visitPoint(const typename Sieve::point_type&) {};
123: inline void visitArrow(const typename Sieve::arrow_type&, const int orientation) {};
124: inline void visitPoint(const typename Sieve::point_type&, const int orientation) {};
125: };
126: class PrintVisitor {
127: protected:
128: ostringstream& os;
129: const int rank;
130: public:
131: PrintVisitor(ostringstream& s, const int rank = 0) : os(s), rank(rank) {};
132: template<typename Arrow>
133: inline void visitArrow(const Arrow& arrow) const {
134: this->os << "["<<this->rank<<"]: " << arrow << std::endl;
135: }
136: template<typename Point>
137: inline void visitPoint(const Point&) const {}
138: };
139: class ReversePrintVisitor : public PrintVisitor {
140: public:
141: ReversePrintVisitor(ostringstream& s, const int rank) : PrintVisitor(s, rank) {};
142: template<typename Arrow>
143: inline void visitArrow(const Arrow& arrow) const {
144: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << std::endl;
145: }
146: template<typename Arrow>
147: inline void visitArrow(const Arrow& arrow, const int orientation) const {
148: this->os << "["<<this->rank<<"]: " << arrow.target << "<----" << arrow.source << ": " << orientation << std::endl;
149: }
150: template<typename Point>
151: inline void visitPoint(const Point&) const {}
152: template<typename Point>
153: inline void visitPoint(const Point&, const int) const {}
154: };
155: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
156: class PointRetriever {
157: public:
158: typedef typename Sieve::point_type point_type;
159: typedef typename Sieve::arrow_type arrow_type;
160: typedef std::pair<point_type,int> oriented_point_type;
161: protected:
162: const bool unique;
163: size_t i, o;
164: size_t skip, limit;
165: Visitor *visitor;
166: size_t size;
167: point_type *points;
168: oriented_point_type *oPoints;
169: protected:
170: inline virtual bool accept(const point_type& point) {return true;};
171: public:
172: PointRetriever() : unique(false), i(0), o(0), skip(0), limit(0) {
173: this->size = 0;
174: this->points = NULL;
175: this->oPoints = NULL;
176: };
177: PointRetriever(const size_t size, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0) {
178: static Visitor nV;
179: this->visitor = &nV;
180: this->points = NULL;
181: this->oPoints = NULL;
182: this->setSize(size);
183: };
184: PointRetriever(const size_t size, Visitor& v, const bool unique = false) : unique(unique), i(0), o(0), skip(0), limit(0), visitor(&v) {
185: this->points = NULL;
186: this->oPoints = NULL;
187: this->setSize(size);
188: };
189: void operator=(const PointRetriever& pr) {
190: this->i = 0; this->o = 0; this->skip = 0; this->limit = 0; this->visitor = pr.visitor;
191: this->points = NULL;
192: this->oPoints = NULL;
193: this->setSize(pr.size);
194: };
195: virtual ~PointRetriever() {
196: delete [] this->points;
197: delete [] this->oPoints;
198: this->points = NULL;
199: this->oPoints = NULL;
200: };
201: inline void visitArrow(const arrow_type& arrow) {
202: this->visitor->visitArrow(arrow);
203: };
204: inline void visitArrow(const arrow_type& arrow, const int orientation) {
205: this->visitor->visitArrow(arrow, orientation);
206: };
207: inline void visitPoint(const point_type& point) {
208: if (i >= size) {
209: ostringstream msg;
210: msg << "Too many points (>" << size << ")for PointRetriever visitor";
211: throw ALE::Exception(msg.str().c_str());
212: }
213: if (this->accept(point)) {
214: if (this->unique) {
215: size_t p;
216: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
217: if (p != i) return;
218: }
219: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
220: points[i++] = point;
221: this->visitor->visitPoint(point);
222: }
223: };
224: inline void visitPoint(const point_type& point, const int orientation) {
225: if (o >= size) {
226: ostringstream msg;
227: msg << "Too many ordered points (>" << size << ")for PointRetriever visitor";
228: throw ALE::Exception(msg.str().c_str());
229: }
230: if (this->accept(point)) {
231: if (this->unique) {
232: size_t p;
233: for(p = 0; p < i; ++p) {if (points[p] == point) break;}
234: if (p != i) return;
235: }
236: if ((i < this->skip) || ((this->limit) && (i >= this->limit))) {--this->skip; return;}
237: points[i++] = point;
238: oPoints[o++] = oriented_point_type(point, orientation);
239: this->visitor->visitPoint(point, orientation);
240: }
241: };
242: public:
243: size_t getSize() const {return this->i;}
244: const point_type *getPoints() const {return this->points;}
245: size_t getOrientedSize() const {return this->o;}
246: const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
247: void clear() {this->i = this->o = 0;}
248: void setSize(const size_t s) {
249: if (this->points) {
250: delete [] this->points;
251: delete [] this->oPoints;
252: }
253: this->size = s;
254: this->points = new point_type[this->size];
255: this->oPoints = new oriented_point_type[this->size];
256: }
257: void setSkip(size_t s) {this->skip = s;};
258: void setLimit(size_t l) {this->limit = l;};
259: };
260: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
261: class NConeRetriever : public PointRetriever<Sieve,Visitor> {
262: public:
263: typedef PointRetriever<Sieve,Visitor> base_type;
264: typedef typename Sieve::point_type point_type;
265: typedef typename Sieve::arrow_type arrow_type;
266: typedef typename base_type::oriented_point_type oriented_point_type;
267: protected:
268: const Sieve& sieve;
269: protected:
270: inline virtual bool accept(const point_type& point) {
271: if (!this->sieve.getConeSize(point))
272: return true;
273: return false;
274: };
275: public:
276: NConeRetriever(const Sieve& s, const size_t size) : PointRetriever<Sieve,Visitor>(size, true), sieve(s) {};
277: NConeRetriever(const Sieve& s, const size_t size, Visitor& v) : PointRetriever<Sieve,Visitor>(size, v, true), sieve(s) {};
278: void operator=(const NConeRetriever& ncr) {
279: this->i = 0; this->o = 0; this->skip = 0; this->limit = 0; this->visitor = ncr.visitor;
280: this->points = NULL;
281: this->oPoints = NULL;
282: this->setSize(ncr.size);
283: };
284: virtual ~NConeRetriever() {};
285: };
286: template<typename Mesh, typename Visitor = NullVisitor<typename Mesh::sieve_type> >
287: class MeshNConeRetriever : public PointRetriever<typename Mesh::sieve_type,Visitor> {
288: public:
289: typedef typename Mesh::Sieve Sieve;
290: typedef PointRetriever<Sieve,Visitor> base_type;
291: typedef typename Sieve::point_type point_type;
292: typedef typename Sieve::arrow_type arrow_type;
293: typedef typename base_type::oriented_point_type oriented_point_type;
294: protected:
295: const Mesh& mesh;
296: const int depth;
297: protected:
298: inline virtual bool accept(const point_type& point) {
299: if (this->mesh.depth(point) == this->depth)
300: return true;
301: return false;
302: };
303: public:
304: MeshNConeRetriever(const Mesh& m, const int depth, const size_t size) : PointRetriever<typename Mesh::Sieve,Visitor>(size, true), mesh(m), depth(depth) {};
305: 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) {};
306: virtual ~MeshNConeRetriever() {};
307: };
308: template<typename Sieve, typename Set, typename Renumbering>
309: class FilteredPointRetriever {
310: public:
311: typedef typename Sieve::point_type point_type;
312: typedef typename Sieve::arrow_type arrow_type;
313: typedef std::pair<point_type,int> oriented_point_type;
314: protected:
315: const Set& pointSet;
316: Renumbering& renumbering;
317: const size_t size;
318: size_t i, o;
319: bool renumber;
320: point_type *points;
321: oriented_point_type *oPoints;
322: public:
323: FilteredPointRetriever(const Set& s, Renumbering& r, const size_t size) : pointSet(s), renumbering(r), size(size), i(0), o(0), renumber(true) {
324: this->points = new point_type[this->size];
325: this->oPoints = new oriented_point_type[this->size];
326: };
327: ~FilteredPointRetriever() {
328: delete [] this->points;
329: delete [] this->oPoints;
330: };
331: inline void visitArrow(const arrow_type& arrow) {};
332: inline void visitPoint(const point_type& point) {
333: if (i >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
334: if (this->pointSet.find(point) == this->pointSet.end()) return;
335: if (renumber) {
336: points[i++] = this->renumbering[point];
337: } else {
338: points[i++] = point;
339: }
340: };
341: inline void visitArrow(const arrow_type& arrow, const int orientation) {};
342: inline void visitPoint(const point_type& point, const int orientation) {
343: if (o >= size) throw ALE::Exception("Too many points for FilteredPointRetriever visitor");
344: if (this->pointSet.find(point) == this->pointSet.end()) return;
345: if (renumber) {
346: points[i++] = this->renumbering[point];
347: oPoints[o++] = oriented_point_type(this->renumbering[point], orientation);
348: } else {
349: points[i++] = point;
350: oPoints[o++] = oriented_point_type(point, orientation);
351: }
352: };
353: public:
354: size_t getSize() const {return this->i;}
355: const point_type *getPoints() const {return this->points;}
356: size_t getOrientedSize() const {return this->o;}
357: const oriented_point_type *getOrientedPoints() const {return this->oPoints;}
358: void clear() {this->i = 0; this->o = 0;}
359: void useRenumbering(const bool renumber) {this->renumber = renumber;}
360: };
361: template<typename Sieve, int size, typename Visitor = NullVisitor<Sieve> >
362: class ArrowRetriever {
363: public:
364: typedef typename Sieve::point_type point_type;
365: typedef typename Sieve::arrow_type arrow_type;
366: typedef std::pair<arrow_type,int> oriented_arrow_type;
367: protected:
368: arrow_type arrows[size];
369: oriented_arrow_type oArrows[size];
370: size_t i, o;
371: Visitor *visitor;
372: public:
373: ArrowRetriever() : i(0), o(0) {
374: static Visitor nV;
375: this->visitor = &nV;
376: };
377: ArrowRetriever(Visitor& v) : i(0), o(0), visitor(&v) {};
378: inline void visitArrow(const typename Sieve::arrow_type& arrow) {
379: if (i >= size) throw ALE::Exception("Too many arrows for visitor");
380: arrows[i++] = arrow;
381: this->visitor->visitArrow(arrow);
382: };
383: inline void visitArrow(const typename Sieve::arrow_type& arrow, const int orientation) {
384: if (o >= size) throw ALE::Exception("Too many arrows for visitor");
385: oArrows[o++] = oriented_arrow_type(arrow, orientation);
386: this->visitor->visitArrow(arrow, orientation);
387: };
388: inline void visitPoint(const point_type& point) {
389: this->visitor->visitPoint(point);
390: };
391: inline void visitPoint(const point_type& point, const int orientation) {
392: this->visitor->visitPoint(point, orientation);
393: };
394: public:
395: size_t getSize() const {return this->i;}
396: const point_type *getArrows() const {return this->arrows;}
397: size_t getOrientedSize() const {return this->o;}
398: const oriented_arrow_type *getOrientedArrows() const {return this->oArrows;}
399: void clear() {this->i = this->o = 0;}
400: };
401: template<typename Sieve, typename Visitor>
402: class ConeVisitor {
403: protected:
404: const Sieve& sieve;
405: Visitor& visitor;
406: bool useSource;
407: public:
408: ConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
409: inline void visitPoint(const typename Sieve::point_type& point) {
410: this->sieve.cone(point, visitor);
411: };
412: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
413: };
414: template<typename Sieve, typename Visitor>
415: class OrientedConeVisitor {
416: protected:
417: const Sieve& sieve;
418: Visitor& visitor;
419: bool useSource;
420: public:
421: OrientedConeVisitor(const Sieve& s, Visitor& v, bool useSource = false) : sieve(s), visitor(v), useSource(useSource) {};
422: inline void visitPoint(const typename Sieve::point_type& point) {
423: this->sieve.orientedCone(point, visitor);
424: };
425: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
426: };
427: template<typename Sieve, typename Visitor>
428: class SupportVisitor {
429: protected:
430: const Sieve& sieve;
431: Visitor& visitor;
432: bool useSource;
433: public:
434: SupportVisitor(const Sieve& s, Visitor& v, bool useSource = true) : sieve(s), visitor(v), useSource(useSource) {};
435: inline void visitPoint(const typename Sieve::point_type& point) {
436: this->sieve.support(point, visitor);
437: };
438: inline void visitArrow(const typename Sieve::arrow_type& arrow) {};
439: };
440: template<typename Sieve, typename Visitor = NullVisitor<Sieve> >
441: class TransitiveClosureVisitor {
442: public:
443: typedef Visitor visitor_type;
444: protected:
445: const Sieve& sieve;
446: Visitor& visitor;
447: bool isCone;
448: std::set<typename Sieve::point_type> seen;
449: public:
450: TransitiveClosureVisitor(const Sieve& s, Visitor& v) : sieve(s), visitor(v), isCone(true) {};
451: inline void visitPoint(const typename Sieve::point_type& point) const {};
452: inline void visitArrow(const typename Sieve::arrow_type& arrow) {
453: if (this->isCone) {
454: if (this->seen.find(arrow.target) == this->seen.end()) {
455: this->seen.insert(arrow.target);
456: this->visitor.visitPoint(arrow.target);
457: }
458: this->visitor.visitArrow(arrow);
459: if (this->seen.find(arrow.source) == this->seen.end()) {
460: if (this->sieve.getConeSize(arrow.source)) {
461: this->sieve.cone(arrow.source, *this);
462: } else {
463: this->seen.insert(arrow.source);
464: this->visitor.visitPoint(arrow.source);
465: }
466: }
467: } else {
468: if (this->seen.find(arrow.source) == this->seen.end()) {
469: this->seen.insert(arrow.source);
470: this->visitor.visitPoint(arrow.source);
471: }
472: this->visitor.visitArrow(arrow);
473: if (this->seen.find(arrow.target) == this->seen.end()) {
474: if (this->sieve.getSupportSize(arrow.target)) {
475: this->sieve.support(arrow.target, *this);
476: } else {
477: this->seen.insert(arrow.target);
478: this->visitor.visitPoint(arrow.target);
479: }
480: }
481: }
482: };
483: public:
484: bool getIsCone() const {return this->isCone;};
485: void setIsCone(const bool isCone) {this->isCone = isCone;};
486: const std::set<typename Sieve::point_type>& getPoints() const {return this->seen;};
487: void clear() {this->seen.clear();};
488: };
489: template<typename Sieve, typename Section>
490: class SizeVisitor {
491: protected:
492: const Section& section;
493: int size;
494: public:
495: SizeVisitor(const Section& s) : section(s), size(0) {};
496: inline void visitPoint(const typename Sieve::point_type& point) {
497: this->size += section.getConstrainedFiberDimension(point);
498: };
499: inline void visitArrow(const typename Sieve::arrow_type&) {};
500: public:
501: int getSize() {return this->size;};
502: };
503: template<typename Sieve, typename Section>
504: class SizeWithBCVisitor {
505: protected:
506: const Section& section;
507: int size;
508: public:
509: SizeWithBCVisitor(const Section& s) : section(s), size(0) {};
510: inline void visitPoint(const typename Sieve::point_type& point) {
511: this->size += section.getFiberDimension(point);
512: };
513: inline void visitArrow(const typename Sieve::arrow_type&) {};
514: public:
515: int getSize() {return this->size;};
516: };
517: template<typename Sieve>
518: class SizeWithBCVisitor<Sieve,PetscSection> {
519: protected:
520: PetscSection section;
521: int size;
522: PetscInt *fieldSize;
523: PetscInt numFields;
524: public:
525: SizeWithBCVisitor(PetscSection s) : section(s), size(0), fieldSize(NULL), numFields(0) {};
526: SizeWithBCVisitor(PetscSection s, PetscInt *fieldSize) : section(s), size(0), fieldSize(fieldSize) {
527: PetscErrorCode PetscSectionGetNumFields(section, &numFields);CHKERRXX(ierr);
528: for(PetscInt f = 0; f < numFields; ++f) {this->fieldSize[f] = 0;}
529: };
530: inline void visitPoint(const typename Sieve::point_type& point) {
531: PetscInt dim;
533: PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
534: this->size += dim;
535: for(PetscInt f = 0; f < numFields; ++f) {
536: PetscSectionGetFieldDof(section, point, f, &dim);CHKERRXX(ierr);
537: this->fieldSize[f] += dim;
538: }
539: };
540: inline void visitArrow(const typename Sieve::arrow_type&) {};
541: public:
542: int getSize() {return this->size;};
543: };
544: template<typename Section>
545: class RestrictVisitor {
546: public:
547: typedef typename Section::value_type value_type;
548: protected:
549: const Section& section;
550: int size;
551: int i;
552: value_type *values;
553: bool allocated;
554: public:
555: RestrictVisitor(const Section& s, const int size) : section(s), size(size), i(0) {
556: this->values = new value_type[this->size];
557: this->allocated = true;
558: };
559: RestrictVisitor(const Section& s, const int size, value_type *values) : section(s), size(size), i(0) {
560: this->values = values;
561: this->allocated = false;
562: };
563: ~RestrictVisitor() {if (this->allocated) {delete [] this->values;}};
564: template<typename Point>
565: inline void visitPoint(const Point& point, const int orientation) {
566: const int dim = section.getFiberDimension(point);
567: if (i+dim > size) {throw ALE::Exception("Too many values for RestrictVisitor.");}
568: const value_type *v = section.restrictPoint(point);
570: if (orientation >= 0) {
571: for(int d = 0; d < dim; ++d, ++i) {
572: this->values[i] = v[d];
573: }
574: } else {
575: for(int d = dim-1; d >= 0; --d, ++i) {
576: this->values[i] = v[d];
577: }
578: }
579: }
580: template<typename Arrow>
581: inline void visitArrow(const Arrow& arrow, const int orientation) {}
582: public:
583: const value_type *getValues() const {return this->values;};
584: int getSize() const {return this->i;};
585: int getMaxSize() const {return this->size;};
586: void ensureSize(const int size) {
587: this->clear();
588: if (size > this->size) {
589: this->size = size;
590: if (this->allocated) {delete [] this->values;}
591: this->values = new value_type[this->size];
592: this->allocated = true;
593: }
594: };
595: void clear() {this->i = 0;};
596: };
597: template<typename ValueType>
598: class RestrictVecVisitor {
599: public:
600: typedef ValueType value_type;
601: protected:
602: const Vec v;
603: const PetscSection section;
604: int size;
605: int i;
606: int nF;
607: int *offsets;
608: int *indices;
609: bool processed;
610: value_type *values;
611: bool allocated;
612: value_type *array;
613: protected:
614: inline void swap(value_type& v, value_type& w) {
615: value_type tmp = v;
616: v = w;
617: w = tmp;
618: };
619: void processArray() {
620: for(PetscInt f = 1; f < nF; ++f) {
621: offsets[f+1] += offsets[f];
622: }
623: for(PetscInt i = 0; i < offsets[nF]; ++i) {
624: indices[i] = offsets[indices[i]]++;
625: }
626: assert(offsets[nF-1] == offsets[nF]);
627: for(PetscInt i = 0; i < offsets[nF]; ++i) {
628: if (indices[i] == -1) continue;
629: PetscInt startPos = indices[i];
630: PetscInt j = startPos, k;
631: value_type val = values[i];
633: do {
634: swap(val, values[j]);
635: k = indices[j];
636: indices[j] = -1;
637: j = k;
638: } while(j != startPos);
639: }
640: };
641: public:
642: RestrictVecVisitor(const Vec v, const PetscSection s, const int size) : v(v), section(s), size(size), i(0), nF(0), processed(true) {
643: this->values = new value_type[this->size];
644: this->allocated = true;
645: PetscErrorCode VecGetArray(this->v, &this->array);CHKERRXX(ierr);
646: };
647: RestrictVecVisitor(const Vec v, const PetscSection s, const int size, value_type *values) : v(v), section(s), size(size), i(0), nF(0), processed(true) {
648: this->values = values;
649: this->allocated = false;
650: PetscErrorCode VecGetArray(this->v, &this->array);CHKERRXX(ierr);
651: };
652: RestrictVecVisitor(const Vec v, const PetscSection s, const int size, value_type *values, int *offsets, int *indices) : v(v), section(s), size(size), i(0), nF(0), offsets(offsets), indices(indices), processed(false) {
655: this->values = values;
656: this->allocated = false;
657: VecGetArray(this->v, &this->array);CHKERRXX(ierr);
658: PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
659: for(PetscInt f = 0; f <= nF; ++f) {offsets[f] = 0;}
660: };
661: ~RestrictVecVisitor() {
662: if (this->allocated) {delete [] this->values;}
663: PetscErrorCode VecRestoreArray(this->v, &this->array);CHKERRXX(ierr);
664: };
665: template<typename Point>
666: inline void visitPoint(const Point& point, const int orientation) {
667: // Known:
668: // Nf: number of fields
669: // Unknown:
670: // Np: number of points
671: // P: points
672: // Algorithm:
673: // Pass 1: Stack up values as before, but also
674: // count size of each field
675: // Comp 1: Sum up sizes to get field offsets
676: // Pass 2: Number each entry with its intended position
677: // Pass 3: Reorder entries
678: // Algorithm if field sizes are known:
679: // Comp 1: Partition array into field components
680: // Pass 1: Stack up values at field offsets
681: PetscInt dim, off;
684: PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
685: if (i+dim > size) {
686: ostringstream msg;
687: msg << "Too many values for RestrictVisitor "<<i+dim<<" > "<<size<< std::endl;
688: throw ALE::Exception(msg.str().c_str());
689: }
690: PetscSectionGetOffset(section, point, &off);CHKERRXX(ierr);
691: const value_type *v = &array[off];
693: if (nF) {
694: for(PetscInt f = 0, fOff = 0; f < nF; ++f) {
695: PetscInt comp, fDim;
697: PetscSectionGetFieldDof(section, point, f, &fDim);CHKERRXX(ierr);
698: offsets[f+1] += fDim;
699: for(PetscInt d = 0; d < fDim; ++d) {
700: indices[i+d] = f;
701: }
702: if (orientation >= 0) {
703: for(PetscInt d = 0; d < fDim; ++d, ++i) {
704: this->values[i] = v[fOff+d];
705: }
706: } else {
707: PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
708: for(PetscInt d = fDim/comp-1; d >= 0; --d) {
709: for(PetscInt c = 0; c < comp; ++c, ++i) {
710: this->values[i] = v[fOff+d*comp+c];
711: }
712: }
713: }
714: fOff += fDim;
715: }
716: } else {
717: if (orientation >= 0) {
718: for(PetscInt d = 0; d < dim; ++d, ++i) {
719: this->values[i] = v[d];
720: }
721: } else {
722: for(PetscInt d = dim-1; d >= 0; --d, ++i) {
723: this->values[i] = v[d];
724: }
725: }
726: }
727: }
728: template<typename Arrow>
729: inline void visitArrow(const Arrow& arrow, const int orientation) {}
730: public:
731: const value_type *getValues() {
732: if (!processed) {processArray(); processed = true;}
733: return this->values;
734: };
735: int getSize() const {return this->i;};
736: int getMaxSize() const {return this->size;};
737: void ensureSize(const int size) {
738: this->clear();
739: if (size > this->size) {
740: this->size = size;
741: if (this->allocated) {delete [] this->values;}
742: this->values = new value_type[this->size];
743: this->allocated = true;
744: }
745: };
746: void clear() {this->i = 0; if (processed) {processed = false;}};
747: };
748: template<typename Section>
749: class UpdateVisitor {
750: public:
751: typedef typename Section::value_type value_type;
752: protected:
753: Section& section;
754: const value_type *values;
755: int i;
756: public:
757: UpdateVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
758: template<typename Point>
759: inline void visitPoint(const Point& point, const int orientation) {
760: const int dim = section.getFiberDimension(point);
761: this->section.updatePoint(point, &this->values[this->i], orientation);
762: this->i += dim;
763: }
764: template<typename Arrow>
765: inline void visitArrow(const Arrow& arrow, const int orientation) {}
766: void clear() {this->i = 0;};
767: };
768: template<typename Section>
769: class UpdateAllVisitor {
770: public:
771: typedef typename Section::value_type value_type;
772: protected:
773: Section& section;
774: const value_type *values;
775: int i;
776: public:
777: UpdateAllVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
778: template<typename Point>
779: inline void visitPoint(const Point& point, const int orientation) {
780: const int dim = section.getFiberDimension(point);
781: this->section.updatePointAll(point, &this->values[this->i], orientation);
782: this->i += dim;
783: }
784: template<typename Arrow>
785: inline void visitArrow(const Arrow& arrow, const int orientation) {}
786: void clear() {this->i = 0;};
787: };
788: template<typename Section>
789: class UpdateAddVisitor {
790: public:
791: typedef typename Section::value_type value_type;
792: protected:
793: Section& section;
794: const value_type *values;
795: int i;
796: public:
797: UpdateAddVisitor(Section& s, const value_type *v) : section(s), values(v), i(0) {};
798: template<typename Point>
799: inline void visitPoint(const Point& point, const int orientation) {
800: const int dim = section.getFiberDimension(point);
801: this->section.updateAddPoint(point, &this->values[this->i], orientation);
802: this->i += dim;
803: }
804: template<typename Arrow>
805: inline void visitArrow(const Arrow& arrow, const int orientation) {}
806: void clear() {this->i = 0;};
807: };
808: template<typename ValueType>
809: class UpdateVecVisitor {
810: public:
811: typedef ValueType value_type;
812: protected:
813: const Vec v;
814: const PetscSection section;
815: const value_type *values;
816: const InsertMode mode;
817: PetscInt nF;
818: PetscInt i;
819: value_type *array;
820: PetscInt *fieldSize;
821: PetscInt *j;
822: protected:
823: inline static void add (value_type& x, value_type y) {x += y;}
824: inline static void insert(value_type& x, value_type y) {x = y;}
825: template<typename Point>
826: void updatePoint(const Point& point, void (*fuse)(value_type&, value_type), const bool setBC, const int orientation = 1) {
827: PetscInt dim; // The number of dof on this point
828: PetscInt cDim; // The nubmer of constraints on this point
829: const PetscInt *cDof; // The indices of the constrained dofs on this point
830: value_type *a; // The values on this point
831: PetscInt offset, cInd = 0;
832: PetscErrorCode ierr;
834: PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
835: PetscSectionGetConstraintDof(section, point, &cDim);CHKERRXX(ierr);
836: PetscSectionGetOffset(section, point, &offset);CHKERRXX(ierr);
837: a = &array[offset];
838: if (!cDim || setBC) {
839: if (orientation >= 0) {
840: for(PetscInt k = 0; k < dim; ++k) {
841: fuse(a[k], values[i+k]);
842: }
843: } else {
844: for(PetscInt k = 0; k < dim; ++k) {
845: fuse(a[k], values[i+dim-k-1]);
846: }
847: }
848: } else {
849: PetscSectionGetConstraintIndices(section, point, &cDof);CHKERRXX(ierr);
850: if (orientation >= 0) {
851: for(PetscInt k = 0; k < dim; ++k) {
852: if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
853: fuse(a[k], values[i+k]);
854: }
855: } else {
856: for(PetscInt k = 0; k < dim; ++k) {
857: if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
858: fuse(a[k], values[i+dim-k-1]);
859: }
860: }
861: }
862: i += dim;
863: }
864: template<typename Point>
865: void updatePointFields(const Point& point, void (*fuse)(value_type&, value_type), const bool setBC, const int orientation = 1) {
866: value_type *a;
867: PetscInt offset;
868: PetscInt fOff = 0;
871: PetscSectionGetOffset(section, point, &offset);CHKERRXX(ierr);
872: a = &array[offset];
873: for(PetscInt f = 0; f < nF; ++f) {
874: PetscInt dim; // The number of dof for field f on this point
875: PetscInt comp; // The number of components for field f on this point
876: PetscInt cDim; // The nubmer of constraints for field f on this point
877: const PetscInt *cDof; // The indices of the constrained dofs for field f on this point
878: PetscInt cInd = 0;
880: PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
881: PetscSectionGetFieldDof(section, point, f, &dim);CHKERRXX(ierr);
882: PetscSectionGetFieldConstraintDof(section, point, f, &cDim);CHKERRXX(ierr);
883: if (!cDim || setBC) {
884: if (orientation >= 0) {
885: for(PetscInt k = 0; k < dim; ++k) {
886: fuse(a[fOff+k], values[j[f]+k]);
887: }
888: } else {
889: for(PetscInt k = dim/comp-1; k >= 0; --k) {
890: for(PetscInt c = 0; c < comp; ++c) {
891: fuse(a[fOff+(dim/comp-1-k)*comp+c], values[j[f]+k*comp+c]);
892: }
893: }
894: }
895: } else {
896: PetscSectionGetFieldConstraintIndices(section, point, f, &cDof);CHKERRXX(ierr);
897: if (orientation >= 0) {
898: for(PetscInt k = 0; k < dim; ++k) {
899: if ((cInd < cDim) && (k == cDof[cInd])) {++cInd; continue;}
900: fuse(a[fOff+k], values[j[f]+k]);
901: }
902: } else {
903: for(PetscInt k = dim/comp-1; k >= 0; --k) {
904: for(PetscInt c = 0; c < comp; ++c) {
905: PetscInt ind = k*comp+c;
906: if ((cInd < cDim) && (ind == cDof[cInd])) {++cInd; continue;}
907: fuse(a[fOff+(dim/comp-1-k)*comp+c], values[j[f]+ind]);
908: }
909: }
910: }
911: }
912: fOff += dim;
913: j[f] += dim;
914: }
915: }
916: public:
917: UpdateVecVisitor(const Vec v, const PetscSection s, const value_type *values, InsertMode mode) : v(v), section(s), values(values), mode(mode), nF(0), i(0) {};
918: UpdateVecVisitor(const Vec v, const PetscSection s, const value_type *values, InsertMode mode, PetscInt numFields, PetscInt fieldSize[]) : v(v), section(s), values(values), mode(mode), nF(numFields), i(0) {
921: VecGetArray(this->v, &this->array);CHKERRXX(ierr);
922: PetscMalloc2(numFields,PetscInt,&this->fieldSize,numFields,PetscInt,&j);CHKERRXX(ierr);
923: for(PetscInt f = 0; f < nF; ++f) {
924: this->fieldSize[f] = fieldSize[f];
925: }
926: this->clear();
927: };
928: ~UpdateVecVisitor() {
930: VecRestoreArray(this->v, &this->array);CHKERRXX(ierr);
931: PetscFree2(fieldSize,j);CHKERRXX(ierr);
932: };
933: template<typename Point>
934: inline void visitPoint(const Point& point, const int orientation) {
935: if (nF) {
936: switch(mode) {
937: case INSERT_VALUES:
938: updatePointFields(point, this->insert, false, orientation);break;
939: case INSERT_ALL_VALUES:
940: updatePointFields(point, this->insert, true, orientation);break;
941: case ADD_VALUES:
942: updatePointFields(point, this->add, false, orientation);break;
943: case ADD_ALL_VALUES:
944: updatePointFields(point, this->add, true, orientation);break;
945: default:
946: throw PETSc::Exception("Invalid mode");
947: }
948: } else {
949: switch(mode) {
950: case INSERT_VALUES:
951: updatePoint(point, this->insert, false, orientation);break;
952: case INSERT_ALL_VALUES:
953: updatePoint(point, this->insert, true, orientation);break;
954: case ADD_VALUES:
955: updatePoint(point, this->add, false, orientation);break;
956: case ADD_ALL_VALUES:
957: updatePoint(point, this->add, true, orientation);break;
958: default:
959: throw PETSc::Exception("Invalid mode");
960: }
961: }
962: }
963: template<typename Arrow>
964: inline void visitArrow(const Arrow& arrow, const int orientation) {}
965: public:
966: void clear() {
967: this->i = 0;
968: if (nF) {
969: j[0] = 0;
970: for(PetscInt f = 1; f < nF; ++f) {
971: j[f] = j[f-1] + fieldSize[f-1];
972: }
973: }
974: };
975: };
976: template<typename Section, typename Order, typename Value>
977: class IndicesVisitor {
978: public:
979: typedef Value value_type;
980: typedef typename Section::point_type point_type;
981: protected:
982: const Section& section;
983: // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
984: Order& order;
985: int size;
986: int i, p;
987: // If false, constrained indices are returned as negative values. Otherwise, they are omitted
988: bool freeOnly;
989: // If true, it allows space for constrained variables (even if the indices are not returned) Wierd
990: bool skipConstraints;
991: value_type *values;
992: bool allocated;
993: point_type *points;
994: bool allocatedPoints;
995: protected:
996: void getUnconstrainedIndices(const point_type& p, const int orientation) {
997: if (i+section.getFiberDimension(p) > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
998: if (orientation >= 0) {
999: const int start = this->order.getIndex(p);
1000: const int end = start + section.getFiberDimension(p);
1002: for(int j = start; j < end; ++j, ++i) {
1003: this->values[i] = j;
1004: }
1005: } else if (!section.getNumSpaces()) {
1006: const int start = this->order.getIndex(p);
1007: const int end = start + section.getFiberDimension(p);
1009: for(int j = end-1; j >= start; --j, ++i) {
1010: this->values[i] = j;
1011: }
1012: } else {
1013: const int numSpaces = section.getNumSpaces();
1014: int start = this->order.getIndex(p);
1016: for(int space = 0; space < numSpaces; ++space) {
1017: const int end = start + section.getFiberDimension(p, space);
1019: for(int j = end-1; j >= start; --j, ++i) {
1020: this->values[i] = j;
1021: }
1022: start = end;
1023: }
1024: }
1025: };
1026: void getConstrainedIndices(const point_type& p, const int orientation) {
1027: const int cDim = this->section.getConstraintDimension(p);
1028: if (i+cDim > size) {throw ALE::Exception("Too many values for IndicesVisitor.");}
1029: typedef typename Section::bc_type::value_type index_type;
1030: const index_type *cDof = this->section.getConstraintDof(p);
1031: const int start = this->order.getIndex(p);
1033: if (orientation >= 0) {
1034: const int dim = this->section.getFiberDimension(p);
1036: for(int j = start, cInd = 0, k = 0; k < dim; ++k) {
1037: if ((cInd < cDim) && (k == cDof[cInd])) {
1038: if (!freeOnly) values[i++] = -(k+1);
1039: if (skipConstraints) ++j;
1040: ++cInd;
1041: } else {
1042: values[i++] = j++;
1043: }
1044: }
1045: } else {
1046: int offset = 0;
1047: int cOffset = 0;
1048: int k = -1;
1050: for(int space = 0; space < section.getNumSpaces(); ++space) {
1051: const int dim = this->section.getFiberDimension(p, space);
1052: const int tDim = this->section.getConstrainedFiberDimension(p, space);
1053: int cInd = (dim - tDim)-1;
1055: k += dim;
1056: for(int l = 0, j = start+tDim+offset; l < dim; ++l, --k) {
1057: if ((cInd >= 0) && (k == cDof[cInd+cOffset])) {
1058: if (!freeOnly) values[i++] = -(offset+l+1);
1059: if (skipConstraints) --j;
1060: --cInd;
1061: } else {
1062: values[i++] = --j;
1063: }
1064: }
1065: k += dim;
1066: offset += dim;
1067: cOffset += dim - tDim;
1068: }
1069: }
1070: };
1071: public:
1072: 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) {
1073: this->values = new value_type[this->size];
1074: this->allocated = true;
1075: if (unique) {
1076: this->points = new point_type[this->size];
1077: this->allocatedPoints = true;
1078: } else {
1079: this->points = NULL;
1080: this->allocatedPoints = false;
1081: }
1082: };
1083: 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) {
1084: this->values = values;
1085: this->allocated = false;
1086: if (unique) {
1087: this->points = new point_type[this->size];
1088: this->allocatedPoints = true;
1089: } else {
1090: this->points = NULL;
1091: this->allocatedPoints = false;
1092: }
1093: };
1094: ~IndicesVisitor() {
1095: if (this->allocated) {delete [] this->values;}
1096: if (this->allocatedPoints) {delete [] this->points;}
1097: };
1098: inline void visitPoint(const point_type& point, const int orientation) {
1099: if (p >= size) {
1100: ostringstream msg;
1101: msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
1102: throw ALE::Exception(msg.str().c_str());
1103: }
1104: if (points) {
1105: int pp;
1106: for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
1107: if (pp != p) return;
1108: points[p++] = point;
1109: }
1110: const int cDim = this->section.getConstraintDimension(point);
1112: if (!cDim) {
1113: this->getUnconstrainedIndices(point, orientation);
1114: } else {
1115: this->getConstrainedIndices(point, orientation);
1116: }
1117: }
1118: template<typename Arrow>
1119: inline void visitArrow(const Arrow& arrow, const int orientation) {}
1120: public:
1121: const value_type *getValues() const {return this->values;};
1122: int getSize() const {return this->i;};
1123: int getMaxSize() const {return this->size;};
1124: void ensureSize(const int size) {
1125: this->clear();
1126: if (size > this->size) {
1127: this->size = size;
1128: if (this->allocated) {delete [] this->values;}
1129: this->values = new value_type[this->size];
1130: this->allocated = true;
1131: if (this->allocatedPoints) {delete [] this->points;}
1132: this->points = new point_type[this->size];
1133: this->allocatedPoints = true;
1134: }
1135: };
1136: void clear() {this->i = 0; this->p = 0;};
1137: };
1138: template<typename Order, typename Value>
1139: class IndicesVisitor<PetscSection, Order, Value> {
1140: public:
1141: typedef Value value_type;
1142: typedef typename Order::point_type point_type;
1143: protected:
1144: const PetscSection& section;
1145: // This can't be const because UniformSection can't have a const restrict(), because of stupid map semantics
1146: Order& order;
1147: int size;
1148: int i, p;
1149: bool setBC; // If true, returns indices for constrained dofs, otherwise negative values are returned
1150: //bool skipConstraints; // If true, do not return constrained indices at all
1151: value_type *values;
1152: bool allocated;
1153: point_type *points;
1154: PetscInt nF;
1155: PetscInt *fieldSize;
1156: PetscInt *j;
1157: protected:
1158: void updatePoint(const point_type& point, const bool setBC, const int orientation = 1) {
1159: PetscInt dim; // The number of dof on this point
1160: PetscInt cDim; // The nubmer of constraints on this point
1161: const PetscInt *cDof; // The indices of the constrained dofs on this point
1162: PetscInt offset = this->order.getIndex(point);
1163: PetscInt cInd = 0;
1164: PetscErrorCode ierr;
1166: PetscSectionGetDof(section, point, &dim);CHKERRXX(ierr);
1167: PetscSectionGetConstraintDof(section, point, &cDim);CHKERRXX(ierr);
1168: if (!cDim || setBC) {
1169: if (orientation >= 0) {
1170: for(PetscInt k = 0; k < dim; ++k) {
1171: values[i+k] = offset+k;
1172: }
1173: } else {
1174: for(PetscInt k = 0; k < dim; ++k) {
1175: values[i+dim-k-1] = offset+k;
1176: }
1177: }
1178: } else {
1179: PetscSectionGetConstraintIndices(section, point, &cDof);CHKERRXX(ierr);
1180: if (orientation >= 0) {
1181: for(PetscInt k = 0; k < dim; ++k) {
1182: if ((cInd < cDim) && (k == cDof[cInd])) {
1183: // Insert check for returning constrained indices
1184: values[i+k] = -(offset+k+1);
1185: ++cInd;
1186: } else {
1187: values[i+k] = offset+k;
1188: }
1189: }
1190: } else {
1191: for(PetscInt k = 0; k < dim; ++k) {
1192: if ((cInd < cDim) && (k == cDof[cInd])) {
1193: // Insert check for returning constrained indices
1194: values[i+dim-k-1] = -(offset+k+1);
1195: ++cInd;
1196: } else {
1197: values[i+dim-k-1] = offset+k;
1198: }
1199: }
1200: }
1201: }
1202: i += dim;
1203: }
1204: void updatePointFields(const point_type& point, const bool setBC, const int orientation = 1) {
1205: PetscInt offset = this->order.getIndex(point);
1206: PetscInt fOff = 0;
1209: for(PetscInt f = 0; f < nF; ++f) {
1210: PetscInt dim; // The number of dof for field f on this point
1211: PetscInt comp; // The number of components for field f on this point
1212: PetscInt cDim; // The nubmer of constraints for field f on this point
1213: const PetscInt *cDof; // The indices of the constrained dofs for field f on this point
1214: PetscInt cInd = 0;
1216: PetscSectionGetFieldComponents(section, f, &comp);CHKERRXX(ierr);
1217: PetscSectionGetFieldDof(section, point, f, &dim);CHKERRXX(ierr);
1218: PetscSectionGetFieldConstraintDof(section, point, f, &cDim);CHKERRXX(ierr);
1219: if (!cDim || setBC) {
1220: if (orientation >= 0) {
1221: for(PetscInt k = 0; k < dim; ++k) {
1222: values[j[f]+k] = offset+fOff+k;
1223: }
1224: } else {
1225: for(PetscInt k = dim/comp-1; k >= 0; --k) {
1226: for(PetscInt c = 0; c < comp; ++c) {
1227: values[j[f]+(dim/comp-1-k)*comp+c] = offset+fOff+k*comp+c;
1228: }
1229: }
1230: }
1231: } else {
1232: PetscSectionGetFieldConstraintIndices(section, point, f, &cDof);CHKERRXX(ierr);
1233: if (orientation >= 0) {
1234: for(PetscInt k = 0; k < dim; ++k) {
1235: if ((cInd < cDim) && (k == cDof[cInd])) {
1236: values[j[f]+k] = -(offset+fOff+k+1);
1237: ++cInd;
1238: } else {
1239: values[j[f]+k] = offset+fOff+k;
1240: }
1241: }
1242: } else {
1243: for(PetscInt k = dim/comp-1; k >= 0; --k) {
1244: for(PetscInt c = 0; c < comp; ++c) {
1245: PetscInt ind = k*comp+c;
1246: if ((cInd < cDim) && (ind == cDof[cInd])) {
1247: values[j[f]+(dim/comp-1-k)*comp+c] = -(offset+fOff+ind+1);
1248: ++cInd;
1249: } else {
1250: values[j[f]+(dim/comp-1-k)*comp+c] = offset+fOff+ind;
1251: }
1252: }
1253: }
1254: }
1255: }
1256: fOff += dim - cDim;
1257: j[f] += dim;
1258: i += dim;
1259: }
1260: }
1261: public:
1262: IndicesVisitor(const PetscSection& s, Order& o, const int size, const bool unique = false, const PetscInt fieldSize[] = NULL) : section(s), order(o), size(size), i(0), p(0), setBC(false) {
1265: PetscMalloc(this->size * sizeof(value_type), &this->values);CHKERRXX(ierr);
1266: this->allocated = true;
1267: this->points = NULL;
1268: if (unique) {
1269: PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
1270: }
1271: nF = 0;
1272: this->fieldSize = this->j = NULL;
1273: if (fieldSize) {
1274: PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
1275: PetscMalloc2(nF,PetscInt,&this->fieldSize,nF,PetscInt,&j);CHKERRXX(ierr);
1276: for(PetscInt f = 0; f < nF; ++f) {
1277: this->fieldSize[f] = fieldSize[f];
1278: }
1279: }
1280: this->clear();
1281: };
1282: IndicesVisitor(const PetscSection& s, Order& o, const int size, value_type *values, const bool unique = false, const PetscInt fieldSize[] = NULL) : section(s), order(o), size(size), i(0), p(0), setBC(false) {
1285: this->values = values;
1286: this->allocated = false;
1287: this->points = NULL;
1288: if (unique) {
1289: PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
1290: }
1291: nF = 0;
1292: this->fieldSize = this->j = NULL;
1293: if (fieldSize) {
1294: PetscSectionGetNumFields(section, &nF);CHKERRXX(ierr);
1295: PetscMalloc2(nF,PetscInt,&fieldSize,nF,PetscInt,&j);CHKERRXX(ierr);
1296: for(PetscInt f = 0; f < nF; ++f) {
1297: this->fieldSize[f] = fieldSize[f];
1298: }
1299: }
1300: this->clear();
1301: };
1302: ~IndicesVisitor() {
1304: if (this->allocated) {PetscFree(values);CHKERRXX(ierr);}
1305: PetscFree(points);CHKERRXX(ierr);
1306: PetscFree2(fieldSize,j);CHKERRXX(ierr);
1307: };
1308: public:
1309: inline void visitPoint(const point_type& point, const int orientation) {
1310: if (p >= size) {
1311: ostringstream msg;
1312: msg << "Too many points (>" << size << ")for IndicesVisitor visitor";
1313: throw ALE::Exception(msg.str().c_str());
1314: }
1315: if (points) {
1316: PetscInt pp;
1317: for(pp = 0; pp < p; ++pp) {if (points[pp] == point) break;}
1318: if (pp != p) return;
1319: points[p++] = point;
1320: }
1321: if (nF) {
1322: updatePointFields(point, setBC, orientation);
1323: } else {
1324: updatePoint(point, setBC, orientation);
1325: }
1326: }
1327: template<typename Arrow>
1328: inline void visitArrow(const Arrow& arrow, const int orientation) {}
1329: public:
1330: const value_type *getValues() const {return this->values;};
1331: int getSize() const {return this->i;};
1332: int getMaxSize() const {return this->size;};
1333: void ensureSize(const int size) {
1334: this->clear();
1335: if (size > this->size) {
1338: this->size = size;
1339: if (this->allocated) {PetscFree(this->values);CHKERRXX(ierr);}
1340: PetscMalloc(this->size * sizeof(value_type), &this->values);CHKERRXX(ierr);
1341: this->allocated = true;
1342: PetscFree(this->points);CHKERRXX(ierr);
1343: PetscMalloc(this->size * sizeof(point_type), &this->points);CHKERRXX(ierr);
1344: }
1345: };
1346: void clear() {
1347: this->p = 0;
1348: this->i = 0;
1349: if (nF) {
1350: j[0] = 0;
1351: for(PetscInt f = 1; f < nF; ++f) {
1352: j[f] = j[f-1] + fieldSize[f-1];
1353: }
1354: }
1355: };
1356: };
1357: template<typename Sieve, typename Label>
1358: class MarkVisitor {
1359: protected:
1360: Label& label;
1361: int marker;
1362: public:
1363: MarkVisitor(Label& l, const int marker) : label(l), marker(marker) {};
1364: inline void visitPoint(const typename Sieve::point_type& point) {
1365: this->label.setCone(this->marker, point);
1366: };
1367: inline void visitArrow(const typename Sieve::arrow_type&) {};
1368: };
1369: };
1371: template<typename Sieve>
1372: class ISieveTraversal {
1373: public:
1374: typedef typename Sieve::point_type point_type;
1375: public:
1376: template<typename Visitor>
1377: static void orientedClosure(const Sieve& sieve, const point_type& p, Visitor& v) {
1378: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
1379: typedef ISieveVisitor::PointRetriever<Sieve,Retriever> TmpRetriever;
1380: Retriever pV(200, v, true); // Correct estimate is pow(std::max(1, sieve->getMaxConeSize()), mesh->depth())
1381: TmpRetriever cV[2] = {TmpRetriever(200,pV), TmpRetriever(200,pV)};
1382: int c = 0;
1384: v.visitPoint(p, 0);
1385: // Cone is guarateed to be ordered correctly
1386: sieve.orientedCone(p, cV[c]);
1388: while(cV[c].getOrientedSize()) {
1389: const typename Retriever::oriented_point_type *cone = cV[c].getOrientedPoints();
1390: const int coneSize = cV[c].getOrientedSize();
1391: c = 1 - c;
1393: for(int p = 0; p < coneSize; ++p) {
1394: const typename Retriever::point_type& point = cone[p].first;
1395: int pO = cone[p].second == 0 ? 1 : cone[p].second;
1396: const int pConeSize = sieve.getConeSize(point);
1398: if (pO < 0) {
1399: if (pO == -pConeSize) {
1400: sieve.orientedReverseCone(point, cV[c]);
1401: } else {
1402: const int numSkip = sieve.getConeSize(point) + pO;
1404: cV[c].setSkip(cV[c].getSize()+numSkip);
1405: cV[c].setLimit(cV[c].getSize()+pConeSize);
1406: sieve.orientedReverseCone(point, cV[c]);
1407: sieve.orientedReverseCone(point, cV[c]);
1408: cV[c].setSkip(0);
1409: cV[c].setLimit(0);
1410: }
1411: } else {
1412: if (pO == 1) {
1413: sieve.orientedCone(point, cV[c]);
1414: } else {
1415: const int numSkip = pO-1;
1417: cV[c].setSkip(cV[c].getSize()+numSkip);
1418: cV[c].setLimit(cV[c].getSize()+pConeSize);
1419: sieve.orientedCone(point, cV[c]);
1420: sieve.orientedCone(point, cV[c]);
1421: cV[c].setSkip(0);
1422: cV[c].setLimit(0);
1423: }
1424: }
1425: }
1426: cV[1-c].clear();
1427: }
1428: #if 0
1429: // These contain arrows paired with orientations from the \emph{previous} arrow
1430: Obj<orientedArrowArray> cone = new orientedArrowArray();
1431: Obj<orientedArrowArray> base = new orientedArrowArray();
1432: coneSet seen;
1434: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
1435: const arrow_type arrow(*c_iter, p);
1437: cone->push_back(oriented_arrow_type(arrow, 1)); // Notice the orientation of leaf cells is always 1
1438: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0])); // However, we return the actual arrow orientation
1439: }
1440: for(int i = 1; i < depth; ++i) {
1441: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
1443: cone->clear();
1444: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
1445: const arrow_type& arrow = b_iter->first; // This is an arrow considered in the previous round
1446: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source); // We are going to get the cone of this guy
1447: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0]; // The orientation of arrow, which is our pO
1449: if (b_iter->second < 0) { // This is the problem. The second orientation is carried up, being from the previous round
1450: o = -(o+1);
1451: }
1452: if (o < 0) {
1453: const int size = pCone->size();
1455: if (o == -size) {
1456: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
1457: if (seen.find(*c_iter) == seen.end()) {
1458: const arrow_type newArrow(*c_iter, arrow.source);
1459: int pointO = orientation->restrictPoint(newArrow)[0];
1461: seen.insert(*c_iter);
1462: cone->push_back(oriented_arrow_type(newArrow, o));
1463: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
1464: }
1465: }
1466: } else {
1467: const int numSkip = size + o;
1468: int count = 0;
1470: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
1471: if (count < numSkip) continue;
1472: if (seen.find(*c_iter) == seen.end()) {
1473: const arrow_type newArrow(*c_iter, arrow.source);
1474: int pointO = orientation->restrictPoint(newArrow)[0];
1476: seen.insert(*c_iter);
1477: cone->push_back(oriented_arrow_type(newArrow, o));
1478: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
1479: }
1480: }
1481: count = 0;
1482: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
1483: if (count >= numSkip) break;
1484: if (seen.find(*c_iter) == seen.end()) {
1485: const arrow_type newArrow(*c_iter, arrow.source);
1486: int pointO = orientation->restrictPoint(newArrow)[0];
1488: seen.insert(*c_iter);
1489: cone->push_back(oriented_arrow_type(newArrow, o));
1490: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
1491: }
1492: }
1493: }
1494: } else {
1495: if (o == 1) {
1496: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
1497: if (seen.find(*c_iter) == seen.end()) {
1498: const arrow_type newArrow(*c_iter, arrow.source);
1500: seen.insert(*c_iter);
1501: cone->push_back(oriented_arrow_type(newArrow, o));
1502: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
1503: }
1504: }
1505: } else {
1506: const int numSkip = o-1;
1507: int count = 0;
1509: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
1510: if (count < numSkip) continue;
1511: if (seen.find(*c_iter) == seen.end()) {
1512: const arrow_type newArrow(*c_iter, arrow.source);
1514: seen.insert(*c_iter);
1515: cone->push_back(oriented_arrow_type(newArrow, o));
1516: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
1517: }
1518: }
1519: count = 0;
1520: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
1521: if (count >= numSkip) break;
1522: if (seen.find(*c_iter) == seen.end()) {
1523: const arrow_type newArrow(*c_iter, arrow.source);
1525: seen.insert(*c_iter);
1526: cone->push_back(oriented_arrow_type(newArrow, o));
1527: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
1528: }
1529: }
1530: }
1531: }
1532: }
1533: }
1534: #endif
1535: }
1536: template<typename Visitor>
1537: static void orientedStar(const Sieve& sieve, const point_type& p, Visitor& v) {
1538: typedef ISieveVisitor::PointRetriever<Sieve,Visitor> Retriever;
1539: Retriever sV[2] = {Retriever(200,v), Retriever(200,v)};
1540: int s = 0;
1542: v.visitPoint(p, 0);
1543: // Support is guarateed to be ordered correctly
1544: sieve.orientedSupport(p, sV[s]);
1546: while(sV[s].getOrientedSize()) {
1547: const typename Retriever::oriented_point_type *support = sV[s].getOrientedPoints();
1548: const int supportSize = sV[s].getOrientedSize();
1549: s = 1 - s;
1551: for(int p = 0; p < supportSize; ++p) {
1552: const typename Retriever::point_type& point = support[p].first;
1553: int pO = support[p].second == 0 ? 1 : support[p].second;
1554: const int pSupportSize = sieve.getSupportSize(point);
1556: if (pO < 0) {
1557: if (pO == -pSupportSize) {
1558: sieve.orientedReverseSupport(point, sV[s]);
1559: } else {
1560: const int numSkip = sieve.getSupportSize(point) + pO;
1562: sV[s].setSkip(sV[s].getSize()+numSkip);
1563: sV[s].setLimit(sV[s].getSize()+pSupportSize);
1564: sieve.orientedReverseSupport(point, sV[s]);
1565: sieve.orientedReverseSupport(point, sV[s]);
1566: sV[s].setSkip(0);
1567: sV[s].setLimit(0);
1568: }
1569: } else {
1570: if (pO == 1) {
1571: sieve.orientedSupport(point, sV[s]);
1572: } else {
1573: const int numSkip = pO-1;
1575: sV[s].setSkip(sV[s].getSize()+numSkip);
1576: sV[s].setLimit(sV[s].getSize()+pSupportSize);
1577: sieve.orientedSupport(point, sV[s]);
1578: sieve.orientedSupport(point, sV[s]);
1579: sV[s].setSkip(0);
1580: sV[s].setLimit(0);
1581: }
1582: }
1583: }
1584: sV[1-s].clear();
1585: }
1586: }
1587: };
1589: namespace IFSieveDef {
1590: template<typename PointType_>
1591: class Sequence {
1592: public:
1593: typedef PointType_ point_type;
1594: class const_iterator {
1595: public:
1596: // Standard iterator typedefs
1597: typedef std::input_iterator_tag iterator_category;
1598: typedef PointType_ value_type;
1599: typedef int difference_type;
1600: typedef value_type* pointer;
1601: typedef value_type& reference;
1602: protected:
1603: const point_type *_data;
1604: int _pos;
1605: public:
1606: const_iterator(const point_type data[], const int pos) : _data(data), _pos(pos) {};
1607: virtual ~const_iterator() {};
1608: public:
1609: virtual bool operator==(const const_iterator& iter) const {return this->_pos == iter._pos;};
1610: virtual bool operator!=(const const_iterator& iter) const {return this->_pos != iter._pos;};
1611: virtual const value_type operator*() const {return this->_data[this->_pos];};
1612: virtual const_iterator& operator++() {++this->_pos; return *this;};
1613: virtual const_iterator operator++(int n) {
1614: const_iterator tmp(*this);
1615: ++this->_pos;
1616: return tmp;
1617: };
1618: };
1619: typedef const_iterator iterator;
1620: protected:
1621: const point_type *_data;
1622: int _begin;
1623: int _end;
1624: public:
1625: Sequence(const point_type data[], const int begin, const int end) : _data(data), _begin(begin), _end(end) {};
1626: virtual ~Sequence() {};
1627: public:
1628: virtual iterator begin() const {return iterator(_data, _begin);};
1629: virtual iterator end() const {return iterator(_data, _end);};
1630: size_t size() const {return(_end - _begin);}
1631: void reset(const point_type data[], const int begin, const int end) {_data = data; _begin = begin; _end = end;};
1632: };
1633: }
1635: // Interval Final Sieve
1636: // This is just two CSR matrices that give cones and supports
1637: // It is completely static and cannot be resized
1638: // It will operator on visitors, rather than sequences (which are messy)
1639: template<typename Point_, typename Allocator_ = malloc_allocator<Point_> >
1640: class IFSieve : public ParallelObject {
1641: public:
1642: // Types
1643: typedef IFSieve<Point_,Allocator_> this_type;
1644: typedef Point_ point_type;
1645: typedef SimpleArrow<point_type,point_type> arrow_type;
1646: typedef typename arrow_type::source_type source_type;
1647: typedef typename arrow_type::target_type target_type;
1648: typedef int index_type;
1649: // Allocators
1650: typedef Allocator_ point_allocator_type;
1651: typedef typename point_allocator_type::template rebind<index_type>::other index_allocator_type;
1652: typedef typename point_allocator_type::template rebind<int>::other int_allocator_type;
1653: // Interval
1654: typedef Interval<point_type, point_allocator_type> chart_type;
1655: // Dynamic structure
1656: typedef std::map<point_type, std::vector<point_type> > newpoints_type;
1657: // Iterator interface
1658: typedef typename IFSieveDef::Sequence<point_type> coneSequence;
1659: typedef typename IFSieveDef::Sequence<point_type> supportSequence;
1660: // Compatibility types for SieveAlgorithms (until we rewrite for visitors)
1661: typedef std::set<point_type> pointSet;
1662: typedef ALE::array<point_type> pointArray;
1663: typedef pointSet coneSet;
1664: typedef pointSet supportSet;
1665: typedef pointArray coneArray;
1666: typedef pointArray supportArray;
1667: protected:
1668: // Arrow Containers
1669: typedef index_type *offsets_type;
1670: typedef point_type *cones_type;
1671: typedef point_type *supports_type;
1672: // Decorators
1673: typedef int *orientations_type;
1674: protected:
1675: // Data
1676: bool indexAllocated;
1677: offsets_type coneOffsets;
1678: offsets_type supportOffsets;
1679: bool pointAllocated;
1680: index_type maxConeSize;
1681: index_type maxSupportSize;
1682: index_type baseSize;
1683: index_type capSize;
1684: cones_type cones;
1685: supports_type supports;
1686: bool orientCones;
1687: orientations_type coneOrientations;
1688: chart_type chart;
1689: int_allocator_type intAlloc;
1690: index_allocator_type indexAlloc;
1691: point_allocator_type pointAlloc;
1692: newpoints_type newCones;
1693: newpoints_type newSupports;
1694: // Sequences
1695: Obj<coneSequence> coneSeq;
1696: Obj<supportSequence> supportSeq;
1697: protected: // Memory Management
1698: void createIndices() {
1699: this->coneOffsets = indexAlloc.allocate(this->chart.size()+1);
1700: this->coneOffsets -= this->chart.min();
1701: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->coneOffsets+i, index_type(0));}
1702: this->supportOffsets = indexAlloc.allocate(this->chart.size()+1);
1703: this->supportOffsets -= this->chart.min();
1704: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(this->supportOffsets+i, index_type(0));}
1705: this->indexAllocated = true;
1706: };
1707: void destroyIndices(const chart_type& chart, offsets_type *coneOffsets, offsets_type *supportOffsets) {
1708: if (*coneOffsets) {
1709: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*coneOffsets)+i);}
1710: *coneOffsets += chart.min();
1711: indexAlloc.deallocate(*coneOffsets, chart.size()+1);
1712: *coneOffsets = NULL;
1713: }
1714: if (*supportOffsets) {
1715: for(index_type i = chart.min(); i <= chart.max(); ++i) {indexAlloc.destroy((*supportOffsets)+i);}
1716: *supportOffsets += chart.min();
1717: indexAlloc.deallocate(*supportOffsets, chart.size()+1);
1718: *supportOffsets = NULL;
1719: }
1720: };
1721: void destroyIndices() {
1722: this->destroyIndices(this->chart, &this->coneOffsets, &this->supportOffsets);
1723: this->indexAllocated = false;
1724: this->maxConeSize = 0;
1725: this->maxSupportSize = 0;
1726: this->baseSize = -1;
1727: this->capSize = -1;
1728: };
1729: void createPoints() {
1730: this->cones = pointAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1731: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->cones+i, point_type(0));}
1732: this->supports = pointAlloc.allocate(this->supportOffsets[this->chart.max()]-this->supportOffsets[this->chart.min()]);
1733: for(index_type i = this->supportOffsets[this->chart.min()]; i < this->supportOffsets[this->chart.max()]; ++i) {pointAlloc.construct(this->supports+i, point_type(0));}
1734: if (orientCones) {
1735: this->coneOrientations = intAlloc.allocate(this->coneOffsets[this->chart.max()]-this->coneOffsets[this->chart.min()]);
1736: for(index_type i = this->coneOffsets[this->chart.min()]; i < this->coneOffsets[this->chart.max()]; ++i) {intAlloc.construct(this->coneOrientations+i, 0);}
1737: }
1738: this->pointAllocated = true;
1739: };
1740: void destroyPoints(const chart_type& chart, const offsets_type coneOffsets, cones_type *cones, const offsets_type supportOffsets, supports_type *supports, orientations_type *coneOrientations) {
1741: if (*cones) {
1742: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*cones)+i);}
1743: pointAlloc.deallocate(*cones, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1744: *cones = NULL;
1745: }
1746: if (*supports) {
1747: for(index_type i = supportOffsets[chart.min()]; i < supportOffsets[chart.max()]; ++i) {pointAlloc.destroy((*supports)+i);}
1748: pointAlloc.deallocate(*supports, supportOffsets[chart.max()]-supportOffsets[chart.min()]);
1749: *supports = NULL;
1750: }
1751: if (*coneOrientations) {
1752: for(index_type i = coneOffsets[chart.min()]; i < coneOffsets[chart.max()]; ++i) {pointAlloc.destroy((*coneOrientations)+i);}
1753: intAlloc.deallocate(*coneOrientations, coneOffsets[chart.max()]-coneOffsets[chart.min()]);
1754: *coneOrientations = NULL;
1755: }
1756: };
1757: void destroyPoints() {
1758: this->destroyPoints(this->chart, this->coneOffsets, &this->cones, this->supportOffsets, &this->supports, &this->coneOrientations);
1759: this->pointAllocated = false;
1760: };
1761: void prefixSum(const offsets_type array) {
1762: for(index_type p = this->chart.min()+1; p <= this->chart.max(); ++p) {
1763: array[p] = array[p] + array[p-1];
1764: }
1765: };
1766: void calculateBaseAndCapSize() {
1767: this->baseSize = 0;
1768: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1769: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1770: ++this->baseSize;
1771: }
1772: }
1773: this->capSize = 0;
1774: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
1775: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
1776: ++this->capSize;
1777: }
1778: }
1779: };
1780: public:
1781: 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) {
1782: this->coneSeq = new coneSequence(NULL, 0, 0);
1783: this->supportSeq = new supportSequence(NULL, 0, 0);
1784: };
1785: 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) {
1786: this->setChart(chart_type(min, max));
1787: this->coneSeq = new coneSequence(NULL, 0, 0);
1788: this->supportSeq = new supportSequence(NULL, 0, 0);
1789: };
1790: ~IFSieve() {
1791: this->destroyPoints();
1792: this->destroyIndices();
1793: };
1794: public: // Accessors
1795: const chart_type& getChart() const {return this->chart;};
1796: void setChart(const chart_type& chart) {
1797: this->destroyPoints();
1798: this->destroyIndices();
1799: this->chart = chart;
1800: this->createIndices();
1801: };
1802: index_type getMaxConeSize() const {
1803: return this->maxConeSize;
1804: };
1805: index_type getMaxSupportSize() const {
1806: return this->maxSupportSize;
1807: };
1808: bool baseContains(const point_type& p) const {
1809: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1810: this->chart.checkPoint(p);
1812: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
1813: return true;
1814: }
1815: return false;
1816: };
1817: bool orientedCones() const {return this->orientCones;};
1818: // Raw array access
1819: offsets_type getConeOffsets() {return this->coneOffsets;};
1820: offsets_type getSupportOffsets() {return this->supportOffsets;};
1821: cones_type getCones() {return this->cones;};
1822: supports_type getSupports() {return this->supports;};
1823: orientations_type getConeOrientations() {return this->coneOrientations;};
1824: public: // Construction
1825: index_type getConeSize(const point_type& p) const {
1826: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1827: this->chart.checkPoint(p);
1828: return this->coneOffsets[p+1]-this->coneOffsets[p];
1829: };
1830: void setConeSize(const point_type& p, const index_type c) {
1831: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1832: this->chart.checkPoint(p);
1833: this->coneOffsets[p+1] = c;
1834: this->maxConeSize = std::max(this->maxConeSize, c);
1835: };
1836: void addConeSize(const point_type& p, const index_type c) {
1837: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1838: this->chart.checkPoint(p);
1839: this->coneOffsets[p+1] += c;
1840: this->maxConeSize = std::max(this->maxConeSize, this->coneOffsets[p+1]);
1841: };
1842: index_type getSupportSize(const point_type& p) const {
1843: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1844: this->chart.checkPoint(p);
1845: return this->supportOffsets[p+1]-this->supportOffsets[p];
1846: };
1847: void setSupportSize(const point_type& p, const index_type s) {
1848: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1849: this->chart.checkPoint(p);
1850: this->supportOffsets[p+1] = s;
1851: this->maxSupportSize = std::max(this->maxSupportSize, s);
1852: };
1853: void addSupportSize(const point_type& p, const index_type s) {
1854: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1855: this->chart.checkPoint(p);
1856: this->supportOffsets[p+1] += s;
1857: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
1858: };
1859: void allocate() {
1860: if (this->pointAllocated) {throw ALE::Exception("IFSieve points have already been allocated.");}
1861: this->prefixSum(this->coneOffsets);
1862: this->prefixSum(this->supportOffsets);
1863: this->createPoints();
1864: this->calculateBaseAndCapSize();
1865: }
1866: // Purely for backwards compatibility
1867: template<typename Color>
1868: void addArrow(const point_type& p, const point_type& q, const Color c, const bool forceSupport = false) {
1869: this->addArrow(p, q, forceSupport);
1870: }
1871: void addArrow(const point_type& p, const point_type& q, const bool forceSupport = false) {
1872: if (!this->chart.hasPoint(q)) {
1873: if (!this->newCones[q].size() && this->chart.hasPoint(q)) {
1874: const index_type start = this->coneOffsets[q];
1875: const index_type end = this->coneOffsets[q+1];
1877: for(int c = start; c < end; ++c) {
1878: this->newCones[q].push_back(this->cones[c]);
1879: }
1880: }
1881: this->newCones[q].push_back(p);
1882: }
1883: if (!this->chart.hasPoint(p) || forceSupport) {
1884: if (!this->newSupports[p].size() && this->chart.hasPoint(p)) {
1885: const index_type start = this->supportOffsets[p];
1886: const index_type end = this->supportOffsets[p+1];
1888: for(int s = start; s < end; ++s) {
1889: this->newSupports[p].push_back(this->supports[s]);
1890: }
1891: }
1892: this->newSupports[p].push_back(q);
1893: }
1894: };
1895: void reallocate() {
1896: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
1897: if (!this->newCones.size() && !this->newSupports.size()) return;
1898: const chart_type oldChart = this->chart;
1899: offsets_type oldConeOffsets = this->coneOffsets;
1900: offsets_type oldSupportOffsets = this->supportOffsets;
1901: cones_type oldCones = this->cones;
1902: supports_type oldSupports = this->supports;
1903: orientations_type oldConeOrientations = this->coneOrientations;
1904: point_type min = this->chart.min();
1905: point_type max = this->chart.max()-1;
1907: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1908: min = std::min(min, c_iter->first);
1909: max = std::max(max, c_iter->first);
1910: }
1911: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1912: min = std::min(min, s_iter->first);
1913: max = std::max(max, s_iter->first);
1914: }
1915: this->chart = chart_type(min, max+1);
1916: this->createIndices();
1917: // Copy sizes (converted from offsets)
1918: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1919: this->coneOffsets[p+1] = oldConeOffsets[p+1]-oldConeOffsets[p];
1920: this->supportOffsets[p+1] = oldSupportOffsets[p+1]-oldSupportOffsets[p];
1921: }
1922: // Inject new sizes
1923: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1924: this->coneOffsets[c_iter->first+1] = c_iter->second.size();
1925: this->maxConeSize = std::max(this->maxConeSize, (int) c_iter->second.size());
1926: }
1927: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1928: this->supportOffsets[s_iter->first+1] = s_iter->second.size();
1929: this->maxSupportSize = std::max(this->maxSupportSize, (int) s_iter->second.size());
1930: }
1931: this->prefixSum(this->coneOffsets);
1932: this->prefixSum(this->supportOffsets);
1933: this->createPoints();
1934: this->calculateBaseAndCapSize();
1935: // Copy cones and supports
1936: for(point_type p = oldChart.min(); p < oldChart.max(); ++p) {
1937: const index_type cStart = this->coneOffsets[p];
1938: const index_type cOStart = oldConeOffsets[p];
1939: const index_type cOEnd = oldConeOffsets[p+1];
1940: const index_type sStart = this->supportOffsets[p];
1941: const index_type sOStart = oldSupportOffsets[p];
1942: const index_type sOEnd = oldSupportOffsets[p+1];
1944: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1945: this->cones[c] = oldCones[cO];
1946: }
1947: for(int sO = sOStart, s = sStart; sO < sOEnd; ++sO, ++s) {
1948: this->supports[s] = oldSupports[sO];
1949: }
1950: if (this->orientCones) {
1951: for(int cO = cOStart, c = cStart; cO < cOEnd; ++cO, ++c) {
1952: this->coneOrientations[c] = oldConeOrientations[cO];
1953: }
1954: }
1955: }
1956: // Inject new cones and supports
1957: for(typename newpoints_type::const_iterator c_iter = this->newCones.begin(); c_iter != this->newCones.end(); ++c_iter) {
1958: index_type start = this->coneOffsets[c_iter->first];
1960: for(typename std::vector<point_type>::const_iterator p_iter = c_iter->second.begin(); p_iter != c_iter->second.end(); ++p_iter) {
1961: this->cones[start++] = *p_iter;
1962: }
1963: if (start != this->coneOffsets[c_iter->first+1]) throw ALE::Exception("Invalid size for new cone array");
1964: }
1965: for(typename newpoints_type::const_iterator s_iter = this->newSupports.begin(); s_iter != this->newSupports.end(); ++s_iter) {
1966: index_type start = this->supportOffsets[s_iter->first];
1968: for(typename std::vector<point_type>::const_iterator p_iter = s_iter->second.begin(); p_iter != s_iter->second.end(); ++p_iter) {
1969: this->supports[start++] = *p_iter;
1970: }
1971: if (start != this->supportOffsets[s_iter->first+1]) throw ALE::Exception("Invalid size for new support array");
1972: }
1973: this->newCones.clear();
1974: this->newSupports.clear();
1975: this->destroyPoints(oldChart, oldConeOffsets, &oldCones, oldSupportOffsets, &oldSupports, &oldConeOrientations);
1976: this->destroyIndices(oldChart, &oldConeOffsets, &oldSupportOffsets);
1977: };
1978: // Recalculate the support offsets and fill the supports
1979: // This is used when an IFSieve is being used as a label
1980: void recalculateLabel() {
1981: ISieveVisitor::PointRetriever<IFSieve> v(1);
1983: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1984: this->supportOffsets[p+1] = 0;
1985: }
1986: this->maxSupportSize = 0;
1987: for(point_type p = this->getChart().min(); p < this->getChart().max(); ++p) {
1988: this->cone(p, v);
1989: if (v.getSize()) {
1990: this->supportOffsets[v.getPoints()[0]+1]++;
1991: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[v.getPoints()[0]+1]);
1992: }
1993: v.clear();
1994: }
1995: this->prefixSum(this->supportOffsets);
1996: this->calculateBaseAndCapSize();
1997: this->symmetrize();
1998: };
1999: void setCone(const point_type cone[], const point_type& p) {
2000: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2001: this->chart.checkPoint(p);
2002: const index_type start = this->coneOffsets[p];
2003: const index_type end = this->coneOffsets[p+1];
2005: for(index_type c = start, i = 0; c < end; ++c, ++i) {
2006: this->cones[c] = cone[i];
2007: }
2008: };
2009: void setCone(const point_type cone, const point_type& p) {
2010: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2011: this->chart.checkPoint(p);
2012: const index_type start = this->coneOffsets[p];
2013: const index_type end = this->coneOffsets[p+1];
2015: if (end - start != 1) {throw ALE::Exception("IFSieve setCone() called with too few points.");}
2016: this->cones[start] = cone;
2017: };
2018: #if 0
2019: template<typename PointSequence>
2020: void setCone(const PointSequence& cone, const point_type& p) {
2021: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2022: this->chart.checkPoint(p);
2023: const index_type start = this->coneOffsets[p];
2024: const index_type end = this->coneOffsets[p+1];
2025: if (cone.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve cone.");}
2026: typename PointSequence::iterator c_iter = cone.begin();
2028: for(index_type c = start; c < end; ++c, ++c_iter) {
2029: this->cones[c] = *c_iter;
2030: }
2031: };
2032: #endif
2033: void setSupport(const point_type& p, const point_type support[]) {
2034: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2035: this->chart.checkPoint(p);
2036: const index_type start = this->supportOffsets[p];
2037: const index_type end = this->supportOffsets[p+1];
2039: for(index_type s = start, i = 0; s < end; ++s, ++i) {
2040: this->supports[s] = support[i];
2041: }
2042: };
2043: #if 0
2044: template<typename PointSequence>
2045: void setSupport(const point_type& p, const PointSequence& support) {
2046: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2047: this->chart.checkPoint(p);
2048: const index_type start = this->supportOffsets[p];
2049: const index_type end = this->supportOffsets[p+1];
2050: if (support.size() != end - start) {throw ALE::Exception("Invalid size for IFSieve support.");}
2051: typename PointSequence::iterator s_iter = support.begin();
2053: for(index_type s = start; s < end; ++s, ++s_iter) {
2054: this->supports[s] = *s_iter;
2055: }
2056: };
2057: #endif
2058: void setConeOrientation(const int coneOrientation[], const point_type& p) {
2059: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2060: this->chart.checkPoint(p);
2061: const index_type start = this->coneOffsets[p];
2062: const index_type end = this->coneOffsets[p+1];
2064: for(index_type c = start, i = 0; c < end; ++c, ++i) {
2065: this->coneOrientations[c] = coneOrientation[i];
2066: }
2067: };
2068: void symmetrizeSizes(const int numCells, const int numCorners, const int cones[], const int offset = 0) {
2069: for(point_type p = 0; p < numCells; ++p) {
2070: const index_type start = p*numCorners;
2071: const index_type end = (p+1)*numCorners;
2073: for(index_type c = start; c < end; ++c) {
2074: const point_type q = cones[c]+offset;
2076: this->supportOffsets[q+1]++;
2077: }
2078: }
2079: for(point_type p = numCells; p < this->chart.max(); ++p) {
2080: this->maxSupportSize = std::max(this->maxSupportSize, this->supportOffsets[p+1]);
2081: }
2082: };
2083: void symmetrize() {
2084: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
2085: offsets -= this->chart.min();
2086: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
2087: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2089: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2090: const index_type start = this->coneOffsets[p];
2091: const index_type end = this->coneOffsets[p+1];
2093: for(index_type c = start; c < end; ++c) {
2094: const point_type q = this->cones[c];
2096: this->chart.checkPoint(q);
2097: this->supports[this->supportOffsets[q]+offsets[q]] = p;
2098: ++offsets[q];
2099: }
2100: }
2101: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
2102: indexAlloc.deallocate(offsets, this->chart.size()+1);
2103: }
2104: index_type getBaseSize() const {
2105: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2106: return this->baseSize;
2107: }
2108: index_type getCapSize() const {
2109: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2110: return this->capSize;
2111: }
2112: template<typename _Section>
2113: void relabel(_Section& labeling) {
2116: index_type *offsets = indexAlloc.allocate(this->chart.size()+1);
2117: offsets -= this->chart.min();
2118: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.construct(offsets+i, index_type(0));}
2119: // Recalculate coneOffsets
2120: for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
2121: const point_type newP = labeling.restrictPoint(p)[0];
2123: offsets[newP+1] = this->getConeSize(p);
2124: }
2125: this->prefixSum(offsets);
2126: PetscMemcpy(this->coneOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
2127: // Recalculate supportOffsets
2128: for(index_type p = this->chart.min(); p < this->chart.max(); ++p) {
2129: const point_type newP = labeling.restrictPoint(p)[0];
2131: offsets[newP+1] = this->getSupportSize(p);
2132: }
2133: this->prefixSum(offsets);
2134: PetscMemcpy(this->supportOffsets, offsets, (this->chart.size()+1)*sizeof(index_type));CHKERRXX(ierr);
2135: for(index_type i = this->chart.min(); i <= this->chart.max(); ++i) {indexAlloc.destroy(offsets+i);}
2136: indexAlloc.deallocate(offsets, this->chart.size()+1);
2137: index_type size = std::max(this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()],
2138: this->supportOffsets[this->chart.max()] - this->supportOffsets[this->chart.min()]);
2139: index_type *orientations = offsets = indexAlloc.allocate(size);
2140: for(index_type i = 0; i < size; ++i) {indexAlloc.construct(orientations+i, index_type(0));}
2141: // Recalculate coneOrientations
2142: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
2143: const point_type newP = labeling.restrictPoint(p)[0];
2144: const index_type start = this->coneOffsets[newP];
2145: const index_type end = this->coneOffsets[newP+1];
2147: for(index_type c = start; c < end; ++c, ++offset) {
2148: orientations[c] = this->coneOrientations[offset];
2149: }
2150: }
2151: PetscMemcpy(this->coneOrientations, orientations, (this->coneOffsets[this->chart.max()] - this->coneOffsets[this->chart.min()])*sizeof(index_type));CHKERRXX(ierr);
2152: for(index_type i = 0; i < size; ++i) {indexAlloc.destroy(orientations+i);}
2153: indexAlloc.deallocate(orientations, size);
2154: // Recalculate cones
2155: point_type *array = pointAlloc.allocate(size);
2157: for(index_type i = 0; i < size; ++i) {pointAlloc.construct(array+i, point_type(0));}
2158: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
2159: const point_type newP = labeling.restrictPoint(p)[0];
2160: const index_type start = this->coneOffsets[newP];
2161: const index_type end = this->coneOffsets[newP+1];
2163: for(index_type c = start; c < end; ++c, ++offset) {
2164: const point_type newQ = labeling.restrictPoint(this->cones[offset])[0];
2166: array[c] = newQ;
2167: }
2168: }
2169: PetscMemcpy(this->cones, array, size*sizeof(point_type));CHKERRXX(ierr);
2170: // Recalculate supports
2171: for(index_type p = this->chart.min(), offset = 0; p < this->chart.max(); ++p) {
2172: const point_type newP = labeling.restrictPoint(p)[0];
2173: const index_type start = this->supportOffsets[newP];
2174: const index_type end = this->supportOffsets[newP+1];
2176: for(index_type c = start; c < end; ++c, ++offset) {
2177: const point_type newQ = labeling.restrictPoint(this->supports[offset])[0];
2179: array[c] = newQ;
2180: }
2181: }
2182: PetscMemcpy(this->supports, array, size*sizeof(point_type));CHKERRXX(ierr);
2183: for(index_type i = 0; i < size; ++i) {pointAlloc.destroy(array+i);}
2184: pointAlloc.deallocate(array, size);
2185: }
2186: public: // Traversals
2187: template<typename Visitor>
2188: void roots(const Visitor& v) const {
2189: this->roots(const_cast<Visitor&>(v));
2190: }
2191: template<typename Visitor>
2192: void roots(Visitor& v) const {
2193: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2195: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2196: if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
2197: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
2198: v.visitPoint(p);
2199: }
2200: }
2201: }
2202: }
2203: int numRoots() {
2204: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2205: int n = 0;
2207: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2208: if (this->coneOffsets[p+1] == this->coneOffsets[p]) {
2209: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
2210: ++n;
2211: }
2212: }
2213: }
2214: return n;
2215: }
2216: template<typename Visitor>
2217: void leaves(const Visitor& v) const {
2218: this->leaves(const_cast<Visitor&>(v));
2219: }
2220: template<typename Visitor>
2221: void leaves(Visitor& v) const {
2222: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2224: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2225: if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
2226: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
2227: v.visitPoint(p);
2228: }
2229: }
2230: }
2231: }
2232: int numLeaves() {
2233: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2234: int n = 0;
2236: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2237: if (this->supportOffsets[p+1] == this->supportOffsets[p]) {
2238: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
2239: ++n;
2240: }
2241: }
2242: }
2243: return n;
2244: }
2245: template<typename Visitor>
2246: void base(const Visitor& v) const {
2247: this->base(const_cast<Visitor&>(v));
2248: }
2249: template<typename Visitor>
2250: void base(Visitor& v) const {
2251: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2253: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2254: if (this->coneOffsets[p+1]-this->coneOffsets[p] > 0) {
2255: v.visitPoint(p);
2256: }
2257: }
2258: }
2259: template<typename Visitor>
2260: void cap(const Visitor& v) const {
2261: this->cap(const_cast<Visitor&>(v));
2262: }
2263: template<typename Visitor>
2264: void cap(Visitor& v) const {
2265: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2267: for(point_type p = this->chart.min(); p < this->chart.max(); ++p) {
2268: if (this->supportOffsets[p+1]-this->supportOffsets[p] > 0) {
2269: v.visitPoint(p);
2270: }
2271: }
2272: }
2273: template<typename PointSequence, typename Visitor>
2274: void cone(const PointSequence& points, Visitor& v) const {
2275: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2276: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2277: const point_type p = *p_iter;
2278: this->chart.checkPoint(p);
2279: const index_type start = this->coneOffsets[p];
2280: const index_type end = this->coneOffsets[p+1];
2282: for(index_type c = start; c < end; ++c) {
2283: v.visitArrow(arrow_type(this->cones[c], p));
2284: v.visitPoint(this->cones[c]);
2285: }
2286: }
2287: }
2288: template<typename Visitor>
2289: void cone(const point_type& p, Visitor& v) const {
2290: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2291: this->chart.checkPoint(p);
2292: const index_type start = this->coneOffsets[p];
2293: const index_type end = this->coneOffsets[p+1];
2295: for(index_type c = start; c < end; ++c) {
2296: v.visitArrow(arrow_type(this->cones[c], p));
2297: v.visitPoint(this->cones[c]);
2298: }
2299: }
2300: const Obj<coneSequence>& cone(const point_type& p) const {
2301: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2302: if (!this->chart.hasPoint(p)) {
2303: this->coneSeq->reset(this->cones, 0, 0);
2304: } else {
2305: this->coneSeq->reset(this->cones, this->coneOffsets[p], this->coneOffsets[p+1]);
2306: }
2307: return this->coneSeq;
2308: }
2309: template<typename PointSequence, typename Visitor>
2310: void support(const PointSequence& points, Visitor& v) const {
2311: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2312: for(typename PointSequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2313: const point_type p = *p_iter;
2314: this->chart.checkPoint(p);
2315: const index_type start = this->supportOffsets[p];
2316: const index_type end = this->supportOffsets[p+1];
2318: for(index_type s = start; s < end; ++s) {
2319: v.visitArrow(arrow_type(p, this->supports[s]));
2320: v.visitPoint(this->supports[s]);
2321: }
2322: }
2323: }
2324: template<typename Visitor>
2325: void support(const point_type& p, Visitor& v) const {
2326: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2327: this->chart.checkPoint(p);
2328: const index_type start = this->supportOffsets[p];
2329: const index_type end = this->supportOffsets[p+1];
2331: for(index_type s = start; s < end; ++s) {
2332: v.visitArrow(arrow_type(p, this->supports[s]));
2333: v.visitPoint(this->supports[s]);
2334: }
2335: }
2336: const Obj<supportSequence>& support(const point_type& p) const {
2337: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2338: if (!this->chart.hasPoint(p)) {
2339: this->supportSeq->reset(this->supports, 0, 0);
2340: } else {
2341: this->supportSeq->reset(this->supports, this->supportOffsets[p], this->supportOffsets[p+1]);
2342: }
2343: return this->supportSeq;
2344: }
2345: template<typename Visitor>
2346: void orientedCone(const point_type& p, Visitor& v) const {
2347: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
2348: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2349: this->chart.checkPoint(p);
2350: const index_type start = this->coneOffsets[p];
2351: const index_type end = this->coneOffsets[p+1];
2353: for(index_type c = start; c < end; ++c) {
2354: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
2355: v.visitPoint(this->cones[c], this->coneOrientations[c]);
2356: }
2357: }
2358: template<typename Visitor>
2359: void orientedReverseCone(const point_type& p, Visitor& v) const {
2360: if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
2361: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2362: this->chart.checkPoint(p);
2363: const index_type start = this->coneOffsets[p];
2364: const index_type end = this->coneOffsets[p+1];
2366: for(index_type c = end-1; c >= start; --c) {
2367: v.visitArrow(arrow_type(this->cones[c], p), this->coneOrientations[c]);
2368: v.visitPoint(this->cones[c], this->coneOrientations[c] ? -(this->coneOrientations[c]+1): 0);
2369: }
2370: }
2371: template<typename Visitor>
2372: void orientedSupport(const point_type& p, Visitor& v) const {
2373: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
2374: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2375: this->chart.checkPoint(p);
2376: const index_type start = this->supportOffsets[p];
2377: const index_type end = this->supportOffsets[p+1];
2379: for(index_type s = start; s < end; ++s) {
2380: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
2381: //v.visitPoint(this->supports[s], this->supportOrientations[s]);
2382: v.visitArrow(arrow_type(this->supports[s], p), 0);
2383: v.visitPoint(this->supports[s], 0);
2384: }
2385: }
2386: template<typename Visitor>
2387: void orientedReverseSupport(const point_type& p, Visitor& v) const {
2388: //if (!this->orientCones) {throw ALE::Exception("IFSieve cones have not been oriented.");}
2389: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2390: this->chart.checkPoint(p);
2391: const index_type start = this->supportOffsets[p];
2392: const index_type end = this->supportOffsets[p+1];
2394: for(index_type s = end-1; s >= start; --s) {
2395: //v.visitArrow(arrow_type(this->supports[s], p), this->supportOrientations[s]);
2396: //v.visitPoint(this->supports[s], this->supportOrientations[s] ? -(this->supportOrientations[s]+1): 0);
2397: v.visitArrow(arrow_type(this->supports[s], p), 0);
2398: v.visitPoint(this->supports[s], 0);
2399: }
2400: }
2401: // Currently does only 1 level
2402: // Does not check for uniqueness
2403: template<typename Visitor>
2404: void meet(const point_type& p, const point_type& q, Visitor& v) const {
2405: if (!this->pointAllocated) {throw ALE::Exception("IFSieve points have not been allocated.");}
2406: this->chart.checkPoint(p);
2407: this->chart.checkPoint(q);
2408: const index_type startP = this->coneOffsets[p];
2409: const index_type endP = this->coneOffsets[p+1];
2410: const index_type startQ = this->coneOffsets[q];
2411: const index_type endQ = this->coneOffsets[q+1];
2413: for(index_type cP = startP; cP < endP; ++cP) {
2414: const point_type& c1 = this->cones[cP];
2416: for(index_type cQ = startQ; cQ < endQ; ++cQ) {
2417: if (c1 == this->cones[cQ]) v.visitPoint(c1);
2418: }
2419: if (this->coneOffsets[c1+1] > this->coneOffsets[c1]) {throw ALE::Exception("Cannot handle multiple level meet()");}
2420: }
2421: }
2422: // Currently does only 1 level
2423: template<typename Sequence, typename Visitor>
2424: void join(const Sequence& points, Visitor& v) const {
2425: typedef std::set<point_type> pointSet;
2426: pointSet intersect[2] = {pointSet(), pointSet()};
2427: pointSet tmp;
2428: int p = 0;
2429: int c = 0;
2431: for(typename Sequence::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
2432: this->chart.checkPoint(*p_iter);
2433: tmp.insert(&this->supports[this->supportOffsets[*p_iter]], &this->supports[this->supportOffsets[(*p_iter)+1]]);
2434: if (p == 0) {
2435: intersect[1-c].insert(tmp.begin(), tmp.end());
2436: p++;
2437: } else {
2438: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
2439: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
2440: intersect[c].clear();
2441: }
2442: c = 1 - c;
2443: tmp.clear();
2444: }
2445: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
2446: v.visitPoint(*p_iter);
2447: }
2448: }
2449: // Helper function
2450: void insertNSupport(point_type p, pointSet& set, const int depth) {
2451: const index_type start = this->supportOffsets[p];
2452: const index_type end = this->supportOffsets[p+1];
2454: if (depth == 1) {
2455: set.insert(&this->supports[start], &this->supports[end]);
2456: } else {
2457: for(index_type s = start; s < end; ++s) {
2458: this->insertNSupport(this->supports[s], set, depth-1);
2459: }
2460: }
2461: }
2462: // Gives only the join of depth n
2463: template<typename SequenceIterator, typename Visitor>
2464: void nJoin(const SequenceIterator& pointsBegin, const SequenceIterator& pointsEnd, const int depth, Visitor& v) {
2465: typedef std::set<point_type> pointSet;
2466: pointSet intersect[2] = {pointSet(), pointSet()};
2467: pointSet tmp;
2468: int p = 0;
2469: int c = 0;
2471: for(SequenceIterator p_iter = pointsBegin; p_iter != pointsEnd; ++p_iter) {
2472: this->chart.checkPoint(*p_iter);
2473: // Put points in the nSupport into tmp (duplicates are fine since it is a set)
2474: this->insertNSupport(*p_iter, tmp, depth);
2475: if (p == 0) {
2476: intersect[1-c].insert(tmp.begin(), tmp.end());
2477: p++;
2478: } else {
2479: std::set_intersection(intersect[c].begin(), intersect[c].end(), tmp.begin(), tmp.end(),
2480: std::insert_iterator<pointSet>(intersect[1-c], intersect[1-c].begin()));
2481: intersect[c].clear();
2482: }
2483: c = 1 - c;
2484: tmp.clear();
2485: }
2486: for(typename pointSet::const_iterator p_iter = intersect[c].begin(); p_iter != intersect[c].end(); ++p_iter) {
2487: v.visitPoint(*p_iter);
2488: }
2489: }
2490: public: // Viewing
2491: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) {
2492: ostringstream txt;
2493: int rank;
2495: if (comm == MPI_COMM_NULL) {
2496: comm = this->comm();
2497: rank = this->commRank();
2498: } else {
2499: MPI_Comm_rank(comm, &rank);
2500: }
2501: if (name == "") {
2502: if(rank == 0) {
2503: txt << "viewing an IFSieve" << std::endl;
2504: }
2505: } else {
2506: if(rank == 0) {
2507: txt << "viewing IFSieve '" << name << "'" << std::endl;
2508: }
2509: }
2510: PetscSynchronizedPrintf(comm, "Max sizes cone: %d support: %d\n", this->getMaxConeSize(), this->getMaxSupportSize());
2511: if(rank == 0) {
2512: txt << "cap --> base:" << std::endl;
2513: }
2514: ISieveVisitor::PrintVisitor pV(txt, rank);
2515: this->cap(ISieveVisitor::SupportVisitor<IFSieve,ISieveVisitor::PrintVisitor>(*this, pV));
2516: PetscSynchronizedPrintf(comm, txt.str().c_str());
2517: PetscSynchronizedFlush(comm);
2518: ostringstream txt2;
2520: if(rank == 0) {
2521: txt2 << "base <-- cap:" << std::endl;
2522: }
2523: ISieveVisitor::ReversePrintVisitor rV(txt2, rank);
2524: this->base(ISieveVisitor::ConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV));
2525: PetscSynchronizedPrintf(comm, txt2.str().c_str());
2526: PetscSynchronizedFlush(comm);
2527: if (orientCones) {
2528: ostringstream txt3;
2530: if(rank == 0) {
2531: txt3 << "Orientation:" << std::endl;
2532: }
2533: ISieveVisitor::ReversePrintVisitor rV2(txt3, rank);
2534: this->base(ISieveVisitor::OrientedConeVisitor<IFSieve,ISieveVisitor::ReversePrintVisitor>(*this, rV2));
2535: PetscSynchronizedPrintf(comm, txt3.str().c_str());
2536: PetscSynchronizedFlush(comm);
2537: }
2538: };
2539: };
2541: class ISieveConverter {
2542: public:
2543: template<typename Sieve, typename ISieve, typename Renumbering>
2544: static void convertSieve(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, bool renumber = true) {
2545: // First construct a renumbering of the sieve points
2546: const Obj<typename Sieve::baseSequence>& base = sieve.base();
2547: const Obj<typename Sieve::capSequence>& cap = sieve.cap();
2548: typename ISieve::point_type min = 0;
2549: typename ISieve::point_type max = 0;
2551: if (renumber) {
2552: /* Roots/Leaves from Sieve do not seem to work */
2554: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2555: if (sieve.support(*b_iter)->size() == 0) {
2556: renumbering[*b_iter] = max++;
2557: }
2558: }
2559: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2560: if (sieve.cone(*c_iter)->size() == 0) {
2561: renumbering[*c_iter] = max++;
2562: }
2563: }
2564: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2565: if (sieve.support(*b_iter)->size() == 0) {
2566: const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
2567: const typename Sieve::coneSequence::iterator coneEnd = sieve.coneEnd(*b_iter);
2569: for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
2570: if (renumbering.find(*c_iter) == renumbering.end()) {
2571: renumbering[*c_iter] = max++;
2572: }
2573: }
2574: }
2575: }
2576: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2577: if (renumbering.find(*b_iter) == renumbering.end()) {
2578: renumbering[*b_iter] = max++;
2579: }
2580: }
2581: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2582: if (renumbering.find(*c_iter) == renumbering.end()) {
2583: renumbering[*c_iter] = max++;
2584: }
2585: }
2586: } else {
2587: if (base->size()) {
2588: min = *base->begin();
2589: max = *base->begin();
2590: } else if (cap->size()) {
2591: min = *cap->begin();
2592: max = *cap->begin();
2593: }
2594: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2595: min = std::min(min, *b_iter);
2596: max = std::max(max, *b_iter);
2597: }
2598: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2599: min = std::min(min, *c_iter);
2600: max = std::max(max, *c_iter);
2601: }
2602: if (base->size() || cap->size()) {
2603: ++max;
2604: }
2605: for(typename ISieve::point_type p = min; p < max; ++p) {
2606: renumbering[p] = p;
2607: }
2608: }
2609: // Create the ISieve
2610: isieve.setChart(typename ISieve::chart_type(min, max));
2611: // Set cone and support sizes
2612: size_t maxSize = 0;
2614: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2615: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2617: isieve.setConeSize(renumbering[*b_iter], cone->size());
2618: maxSize = std::max(maxSize, cone->size());
2619: }
2620: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2621: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2623: isieve.setSupportSize(renumbering[*c_iter], support->size());
2624: maxSize = std::max(maxSize, support->size());
2625: }
2626: isieve.allocate();
2627: // Fill up cones and supports
2628: typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];
2630: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2631: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2632: int i = 0;
2634: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2635: points[i] = renumbering[*c_iter];
2636: }
2637: isieve.setCone(points, renumbering[*b_iter]);
2638: }
2639: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2640: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2641: int i = 0;
2643: for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
2644: points[i] = renumbering[*s_iter];
2645: }
2646: isieve.setSupport(renumbering[*c_iter], points);
2647: }
2648: delete [] points;
2649: }
2650: template<typename Sieve, typename Renumbering>
2651: static void convertSieve(Sieve& sieve, DM dm, Renumbering& renumbering, bool renumber = true) {
2652: // First construct a renumbering of the sieve points
2653: const Obj<typename Sieve::baseSequence>& base = sieve.base();
2654: const Obj<typename Sieve::capSequence>& cap = sieve.cap();
2655: PetscInt min = 0;
2656: PetscInt max = 0;
2657: PetscErrorCode ierr;
2659: if (renumber) {
2660: /* Roots/Leaves from Sieve do not seem to work */
2662: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2663: if (sieve.support(*b_iter)->size() == 0) {
2664: renumbering[*b_iter] = max++;
2665: }
2666: }
2667: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2668: if (sieve.cone(*c_iter)->size() == 0) {
2669: renumbering[*c_iter] = max++;
2670: }
2671: }
2672: #if 0
2673: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2674: if (sieve.support(*b_iter)->size() == 0) {
2675: const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
2676: const typename Sieve::coneSequence::iterator coneEnd = sieve.coneEnd(*b_iter);
2678: for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
2679: if (renumbering.find(*c_iter) == renumbering.end()) {
2680: renumbering[*c_iter] = max++;
2681: }
2682: }
2683: }
2684: }
2685: #else
2686: std::vector<typename Sieve::point_type> faces;
2687: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2688: if (sieve.support(*b_iter)->size() == 0) {
2689: const typename Sieve::coneSequence::iterator coneBegin = sieve.coneBegin(*b_iter);
2690: const typename Sieve::coneSequence::iterator coneEnd = sieve.coneEnd(*b_iter);
2692: for(typename Sieve::coneSequence::iterator c_iter = coneBegin; c_iter != coneEnd; ++c_iter) {
2693: if (renumbering.find(*c_iter) == renumbering.end()) {
2694: faces.push_back(*c_iter);
2695: }
2696: }
2697: }
2698: }
2699: std::sort(faces.begin(), faces.end());
2700: typename std::vector<typename Sieve::point_type>::const_iterator fEnd = std::unique(faces.begin(), faces.end());
2701: for(typename std::vector<typename Sieve::point_type>::const_iterator c_iter = faces.begin(); c_iter != fEnd; ++c_iter) {
2702: renumbering[*c_iter] = max++;
2703: }
2704: faces.clear();
2705: #endif
2706: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2707: if (renumbering.find(*b_iter) == renumbering.end()) {
2708: renumbering[*b_iter] = max++;
2709: }
2710: }
2711: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2712: if (renumbering.find(*c_iter) == renumbering.end()) {
2713: renumbering[*c_iter] = max++;
2714: }
2715: }
2716: } else {
2717: if (base->size()) {
2718: min = *base->begin();
2719: max = *base->begin();
2720: } else if (cap->size()) {
2721: min = *cap->begin();
2722: max = *cap->begin();
2723: }
2724: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2725: min = std::min(min, *b_iter);
2726: max = std::max(max, *b_iter);
2727: }
2728: for(typename Sieve::baseSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2729: min = std::min(min, *c_iter);
2730: max = std::max(max, *c_iter);
2731: }
2732: if (base->size() || cap->size()) {
2733: ++max;
2734: }
2735: for(PetscInt p = min; p < max; ++p) {
2736: renumbering[p] = p;
2737: }
2738: }
2739: // Create the ISieve
2740: DMPlexSetChart(dm, min, max);CHKERRXX(ierr);
2741: // Set cone and support sizes
2742: size_t maxSize = 0;
2744: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2745: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2747: DMPlexSetConeSize(dm, renumbering[*b_iter], cone->size());CHKERRXX(ierr);
2748: maxSize = std::max(maxSize, cone->size());
2749: }
2750: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2751: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2753: DMPlexSetSupportSize(dm, renumbering[*c_iter], support->size());CHKERRXX(ierr);
2754: maxSize = std::max(maxSize, support->size());
2755: }
2756: DMSetUp(dm);CHKERRXX(ierr);
2757: // Fill up cones and supports
2758: typename Sieve::point_type *points = new typename Sieve::point_type[maxSize];
2760: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2761: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2762: int i = 0;
2764: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2765: points[i] = renumbering[*c_iter];
2766: }
2767: DMPlexSetCone(dm, renumbering[*b_iter], points);CHKERRXX(ierr);
2768: }
2769: for(typename Sieve::capSequence::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
2770: const Obj<typename Sieve::supportSequence>& support = sieve.support(*c_iter);
2771: int i = 0;
2773: for(typename Sieve::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter, ++i) {
2774: points[i] = renumbering[*s_iter];
2775: }
2776: DMPlexSetSupport(dm, renumbering[*c_iter], points);CHKERRXX(ierr);
2777: }
2778: delete [] points;
2779: }
2780: template<typename Sieve, typename ISieve, typename Renumbering, typename ArrowSection>
2781: static void convertOrientation(Sieve& sieve, ISieve& isieve, Renumbering& renumbering, ArrowSection *orientation) {
2782: if (isieve.getMaxConeSize() < 0) return;
2783: const Obj<typename Sieve::baseSequence>& base = sieve.base();
2784: int *orientations = new int[isieve.getMaxConeSize()];
2786: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2787: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2788: int i = 0;
2790: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2791: typename ArrowSection::point_type arrow(*c_iter, *b_iter);
2793: orientations[i] = orientation->restrictPoint(arrow)[0];
2794: }
2795: isieve.setConeOrientation(orientations, renumbering[*b_iter]);
2796: }
2797: delete [] orientations;
2798: }
2799: template<typename Sieve, typename Renumbering, typename ArrowSection>
2800: static void convertOrientation(Sieve& sieve, DM dm, Renumbering& renumbering, ArrowSection *orientation) {
2801: PetscInt maxConeSize;
2804: DMPlexGetMaxSizes(dm, &maxConeSize, NULL);CHKERRXX(ierr);
2805: if (maxConeSize < 0) return;
2806: const Obj<typename Sieve::baseSequence>& base = sieve.base();
2807: int *orientations = new int[maxConeSize];
2809: for(typename Sieve::baseSequence::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
2810: const Obj<typename Sieve::coneSequence>& cone = sieve.cone(*b_iter);
2811: int i = 0;
2813: for(typename Sieve::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter, ++i) {
2814: typename ArrowSection::point_type arrow(*c_iter, *b_iter);
2815: const int o = orientation->restrictPoint(arrow)[0];
2817: orientations[i] = o == 1 ? 0 : o;
2818: }
2819: DMPlexSetConeOrientation(dm, renumbering[*b_iter], orientations);
2820: }
2821: delete [] orientations;
2822: }
2823: template<typename Section, typename ISection, typename Renumbering>
2824: static void convertCoordinates(Section& coordinates, ISection& icoordinates, Renumbering& renumbering) {
2825: const typename Section::chart_type& chart = coordinates.getChart();
2826: typename ISection::point_type min = *chart.begin();
2827: typename ISection::point_type max = *chart.begin();
2829: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2830: min = std::min(min, renumbering[*p_iter]);
2831: max = std::max(max, renumbering[*p_iter]);
2832: }
2833: icoordinates.setChart(typename ISection::chart_type(min, max+1));
2834: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2835: icoordinates.setFiberDimension(renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));
2836: }
2837: icoordinates.allocatePoint();
2838: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2839: icoordinates.updatePoint(renumbering[*p_iter], coordinates.restrictPoint(*p_iter));
2840: }
2841: }
2842: template<typename Section, typename Renumbering>
2843: static void convertCoordinates(Section& coordinates, PetscSection coordSection, Vec coords, Renumbering& renumbering) {
2844: const typename Section::chart_type& chart = coordinates.getChart();
2845: PetscInt min = *chart.begin();
2846: PetscInt max = *chart.begin();
2847: PetscScalar *a;
2848: PetscInt n;
2849: PetscErrorCode ierr;
2851: PetscSectionSetNumFields(coordSection, 1);CHKERRXX(ierr);
2852: if (!chart.size()) {
2853: PetscSectionSetFieldComponents(coordSection, 0, 1);CHKERRXX(ierr);
2854: return;
2855: }
2856: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2857: min = std::min(min, renumbering[*p_iter]);
2858: max = std::max(max, renumbering[*p_iter]);
2859: }
2860: PetscSectionSetFieldComponents(coordSection, 0, coordinates.getFiberDimension(*chart.begin()));CHKERRXX(ierr);
2861: PetscSectionSetChart(coordSection, min, max+1);CHKERRXX(ierr);
2862: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2863: PetscSectionSetDof(coordSection, renumbering[*p_iter], coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
2864: PetscSectionSetFieldDof(coordSection, renumbering[*p_iter], 0, coordinates.getFiberDimension(*p_iter));CHKERRXX(ierr);
2865: }
2866: PetscSectionSetUp(coordSection);CHKERRXX(ierr);
2867: PetscSectionGetStorageSize(coordSection, &n);CHKERRXX(ierr);
2868: VecSetSizes(coords, n, PETSC_DETERMINE);CHKERRXX(ierr);
2869: VecSetFromOptions(coords);CHKERRXX(ierr);
2870: VecGetArray(coords, &a);CHKERRXX(ierr);
2871: for(typename Section::chart_type::const_iterator p_iter = chart.begin(); p_iter != chart.end(); ++p_iter) {
2872: const typename Section::value_type *values = coordinates.restrictPoint(*p_iter);
2873: PetscInt dof, off;
2875: PetscSectionGetDof(coordSection, renumbering[*p_iter], &dof);CHKERRXX(ierr);
2876: PetscSectionGetOffset(coordSection, renumbering[*p_iter], &off);CHKERRXX(ierr);
2877: for(int d = 0; d < dof; ++d) {
2878: a[off+d] = values[d];
2879: }
2880: }
2881: VecRestoreArray(coords, &a);CHKERRXX(ierr);
2882: }
2883: template<typename Label, typename Renumbering>
2884: static void convertLabel(const Obj<Label>& newLabel, const Obj<Label>& oldLabel, Renumbering& renumbering) {
2885: for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
2886: if (oldLabel->getConeSize(p->first)) {
2887: newLabel->setCone(*oldLabel->cone(p->first)->begin(), p->second);
2888: }
2889: }
2890: }
2891: template<typename Label, typename Renumbering>
2892: static void convertLabel(DM dm, const char name[], const Obj<Label>& label, Renumbering& renumbering) {
2895: for(typename Renumbering::const_iterator p = renumbering.begin(); p != renumbering.end(); ++p) {
2896: if (label->getConeSize(p->first)) {
2897: DMPlexSetLabelValue(dm, name, p->second, *label->cone(p->first)->begin());CHKERRXX(ierr);
2898: }
2899: }
2900: }
2901: template<typename Mesh, typename IMesh, typename Renumbering>
2902: static void convertMesh(Mesh& mesh, IMesh& imesh, Renumbering& renumbering, bool renumber = true) {
2903: convertSieve(*mesh.getSieve(), *imesh.getSieve(), renumbering, renumber);
2904: imesh.stratify();
2905: convertOrientation(*mesh.getSieve(), *imesh.getSieve(), renumbering, mesh.getArrowSection("orientation").ptr());
2906: convertCoordinates(*mesh.getRealSection("coordinates"), *imesh.getRealSection("coordinates"), renumbering);
2907: if (mesh.hasRealSection("normals")) {
2908: convertCoordinates(*mesh.getRealSection("normals"), *imesh.getRealSection("normals"), renumbering);
2909: }
2910: const typename Mesh::labels_type& labels = mesh.getLabels();
2911: std::string heightName("height");
2912: std::string depthName("depth");
2914: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
2915: #ifdef IMESH_NEW_LABELS
2916: if ((l_iter->first != heightName) && (l_iter->first != depthName)) {
2917: convertLabel(imesh, l_iter->first, l_iter->second);
2918: }
2919: #else
2920: if (renumber) {
2921: convertLabel(imesh.createLabel(l_iter->first), l_iter->second, renumbering);
2922: } else {
2923: imesh.setLabel(l_iter->first, l_iter->second);
2924: }
2925: #endif
2926: }
2927: }
2928: template<typename Mesh, typename Renumbering>
2929: static void convertMesh(Mesh& mesh, DM *dm, Renumbering& renumbering, bool renumber = true) {
2930: PetscSection coordSection;
2931: Vec coordinates;
2934: DMCreate(mesh.comm(), dm);CHKERRXX(ierr);
2935: DMSetType(*dm, DMPLEX);CHKERRXX(ierr);
2936: DMPlexSetDimension(*dm, mesh.getDimension());CHKERRXX(ierr);
2937: convertSieve(*mesh.getSieve(), *dm, renumbering, renumber);
2938: DMPlexStratify(*dm);CHKERRXX(ierr);
2939: convertOrientation(*mesh.getSieve(), *dm, renumbering, mesh.getArrowSection("orientation").ptr());
2940: DMPlexGetCoordinateSection(*dm, &coordSection);CHKERRXX(ierr);
2941: VecCreate(mesh.comm(), &coordinates);CHKERRXX(ierr);
2942: convertCoordinates(*mesh.getRealSection("coordinates"), coordSection, coordinates, renumbering);
2943: DMSetCoordinatesLocal(*dm, coordinates);CHKERRXX(ierr);
2944: VecDestroy(&coordinates);CHKERRXX(ierr);
2945: const typename Mesh::labels_type& labels = mesh.getLabels();
2947: for(typename Mesh::labels_type::const_iterator l_iter = labels.begin(); l_iter != labels.end(); ++l_iter) {
2948: convertLabel(*dm, l_iter->first.c_str(), l_iter->second, renumbering);
2949: }
2950: }
2951: };
2953: class ISieveSerializer {
2954: public:
2955: template<typename ISieve>
2956: static void writeSieve(const std::string& filename, ISieve& sieve) {
2957: std::ofstream fs;
2959: if (sieve.commRank() == 0) {
2960: fs.open(filename.c_str());
2961: }
2962: writeSieve(fs, sieve);
2963: if (sieve.commRank() == 0) {
2964: fs.close();
2965: }
2966: }
2967: template<typename ISieve>
2968: static void writeSieve(std::ofstream& fs, ISieve& sieve) {
2969: typedef ISieveVisitor::PointRetriever<ISieve> Visitor;
2970: const Obj<typename ISieve::chart_type>& chart = sieve.getChart();
2971: typename ISieve::point_type min = chart->min();
2972: typename ISieve::point_type max = chart->max();
2973: PetscInt *mins = new PetscInt[sieve.commSize()];
2974: PetscInt *maxs = new PetscInt[sieve.commSize()];
2976: // Write sizes
2977: if (sieve.commRank() == 0) {
2978: // Write local
2979: fs << min <<" "<< max << std::endl;
2980: for(typename ISieve::point_type p = min; p < max; ++p) {
2981: fs << sieve.getConeSize(p) << " " << sieve.getSupportSize(p) << std::endl;
2982: }
2983: // Receive and write remote
2984: for(int pr = 1; pr < sieve.commSize(); ++pr) {
2985: PetscInt minp, maxp;
2986: PetscInt *sizes;
2987: MPI_Status status;
2990: MPI_Recv(&minp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2991: MPI_Recv(&maxp, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2992: PetscMalloc(2*(maxp - minp) * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
2993: MPI_Recv(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
2994: fs << minp <<" "<< maxp << std::endl;
2995: for(PetscInt s = 0; s < maxp - minp; ++s) {
2996: fs << sizes[s*2+0] << " " << sizes[s*2+1] << std::endl;
2997: }
2998: PetscFree(sizes);CHKERRXX(ierr);
2999: mins[pr] = minp;
3000: maxs[pr] = maxp;
3001: }
3002: } else {
3003: // Send remote
3004: //PetscInt min = chart->min();
3005: //PetscInt max = chart->max();
3006: PetscInt s = 0;
3007: PetscInt *sizes;
3010: MPI_Send(&min, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
3011: MPI_Send(&max, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
3012: PetscMalloc((max - min)*2 * sizeof(PetscInt) + 1, &sizes);CHKERRXX(ierr);
3013: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
3014: sizes[s*2+0] = sieve.getConeSize(p);
3015: sizes[s*2+1] = sieve.getSupportSize(p);
3016: }
3017: MPI_Send(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
3018: PetscFree(sizes);CHKERRXX(ierr);
3019: }
3020: // Write covering
3021: if (sieve.commRank() == 0) {
3022: // Write local
3023: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
3025: for(typename ISieve::point_type p = min; p < max; ++p) {
3026: sieve.cone(p, pV);
3027: const typename Visitor::point_type *cone = pV.getPoints();
3028: const int cSize = pV.getSize();
3030: if (cSize > 0) {
3031: for(int c = 0; c < cSize; ++c) {
3032: if (c) {fs << " ";}
3033: fs << cone[c];
3034: }
3035: fs << std::endl;
3036: }
3037: pV.clear();
3039: sieve.orientedCone(p, pV);
3040: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
3041: const int oSize = pV.getOrientedSize();
3043: if (oSize > 0) {
3044: for(int c = 0; c < oSize; ++c) {
3045: if (c) {fs << " ";}
3046: fs << oCone[c].second;
3047: }
3048: fs << std::endl;
3049: }
3050: pV.clear();
3052: sieve.support(p, pV);
3053: const typename Visitor::point_type *support = pV.getPoints();
3054: const int sSize = pV.getSize();
3056: if (sSize > 0) {
3057: for(int s = 0; s < sSize; ++s) {
3058: if (s) {fs << " ";}
3059: fs << support[s];
3060: }
3061: fs << std::endl;
3062: }
3063: pV.clear();
3064: }
3065: // Receive and write remote
3066: for(int pr = 1; pr < sieve.commSize(); ++pr) {
3067: PetscInt size;
3068: PetscInt *data;
3069: PetscInt off = 0;
3070: MPI_Status status;
3073: MPI_Recv(&size, 1, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
3074: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
3075: MPI_Recv(data, size, MPIU_INT, pr, 1, sieve.comm(), &status);CHKERRXX(ierr);
3076: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
3077: PetscInt cSize = data[off++];
3079: fs << cSize << std::endl;
3080: if (cSize > 0) {
3081: for(int c = 0; c < cSize; ++c) {
3082: if (c) {fs << " ";}
3083: fs << data[off++];
3084: }
3085: fs << std::endl;
3086: }
3087: PetscInt oSize = data[off++];
3089: if (oSize > 0) {
3090: for(int c = 0; c < oSize; ++c) {
3091: if (c) {fs << " ";}
3092: fs << data[off++];
3093: }
3094: fs << std::endl;
3095: }
3096: PetscInt sSize = data[off++];
3098: fs << sSize << std::endl;
3099: if (sSize > 0) {
3100: for(int s = 0; s < sSize; ++s) {
3101: if (s) {fs << " ";}
3102: fs << data[off++];
3103: }
3104: fs << std::endl;
3105: }
3106: }
3107: assert(off == size);
3108: PetscFree(data);CHKERRXX(ierr);
3109: }
3110: } else {
3111: // Send remote
3112: Visitor pV(std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize())+1);
3113: PetscInt totalConeSize = 0;
3114: PetscInt totalSupportSize = 0;
3116: for(typename ISieve::point_type p = min; p < max; ++p) {
3117: totalConeSize += sieve.getConeSize(p);
3118: totalSupportSize += sieve.getSupportSize(p);
3119: }
3120: PetscInt size = (sieve.getChart().size()+totalConeSize)*2 + sieve.getChart().size()+totalSupportSize;
3121: PetscInt off = 0;
3122: PetscInt *data;
3125: MPI_Send(&size, 1, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
3126: // There is no nice way to make a generic MPI type here. Really sucky
3127: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
3128: for(typename ISieve::point_type p = min; p < max; ++p) {
3129: sieve.cone(p, pV);
3130: const typename Visitor::point_type *cone = pV.getPoints();
3131: const int cSize = pV.getSize();
3133: data[off++] = cSize;
3134: for(int c = 0; c < cSize; ++c) {
3135: data[off++] = cone[c];
3136: }
3137: pV.clear();
3139: sieve.orientedCone(p, pV);
3140: const typename Visitor::oriented_point_type *oCone = pV.getOrientedPoints();
3141: const int oSize = pV.getOrientedSize();
3143: data[off++] = oSize;
3144: for(int c = 0; c < oSize; ++c) {
3145: data[off++] = oCone[c].second;
3146: }
3147: pV.clear();
3149: sieve.support(p, pV);
3150: const typename Visitor::point_type *support = pV.getPoints();
3151: const int sSize = pV.getSize();
3153: data[off++] = sSize;
3154: for(int s = 0; s < sSize; ++s) {
3155: data[off++] = support[s];
3156: }
3157: pV.clear();
3158: }
3159: MPI_Send(data, size, MPIU_INT, 0, 1, sieve.comm());CHKERRXX(ierr);
3160: PetscFree(data);CHKERRXX(ierr);
3161: }
3162: delete [] mins;
3163: delete [] maxs;
3164: // Output renumbering
3165: }
3166: template<typename ISieve>
3167: static void loadSieve(const std::string& filename, ISieve& sieve) {
3168: std::ifstream fs;
3170: if (sieve.commRank() == 0) {
3171: fs.open(filename.c_str());
3172: }
3173: loadSieve(fs, sieve);
3174: if (sieve.commRank() == 0) {
3175: fs.close();
3176: }
3177: }
3178: template<typename ISieve>
3179: static void loadSieve(std::ifstream& fs, ISieve& sieve) {
3180: typename ISieve::point_type min, max;
3181: PetscInt *mins = new PetscInt[sieve.commSize()];
3182: PetscInt *maxs = new PetscInt[sieve.commSize()];
3183: PetscInt *totalConeSizes = new PetscInt[sieve.commSize()];
3184: PetscInt *totalSupportSizes = new PetscInt[sieve.commSize()];
3186: // Load sizes
3187: if (sieve.commRank() == 0) {
3188: // Load local
3189: fs >> min;
3190: fs >> max;
3191: sieve.setChart(typename ISieve::chart_type(min, max));
3192: for(typename ISieve::point_type p = min; p < max; ++p) {
3193: typename ISieve::index_type coneSize, supportSize;
3195: fs >> coneSize;
3196: fs >> supportSize;
3197: sieve.setConeSize(p, coneSize);
3198: sieve.setSupportSize(p, supportSize);
3199: }
3200: // Load and send remote
3201: for(int pr = 1; pr < sieve.commSize(); ++pr) {
3202: PetscInt minp, maxp;
3203: PetscInt *sizes;
3206: fs >> minp;
3207: fs >> maxp;
3208: MPI_Send(&minp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
3209: MPI_Send(&maxp, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
3210: PetscMalloc((maxp - minp)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
3211: totalConeSizes[pr] = 0;
3212: totalSupportSizes[pr] = 0;
3213: for(PetscInt s = 0; s < maxp - minp; ++s) {
3214: fs >> sizes[s*2+0];
3215: fs >> sizes[s*2+1];
3216: totalConeSizes[pr] += sizes[s*2+0];
3217: totalSupportSizes[pr] += sizes[s*2+1];
3218: }
3219: MPI_Send(sizes, (maxp - minp)*2, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
3220: PetscFree(sizes);CHKERRXX(ierr);
3221: mins[pr] = minp;
3222: maxs[pr] = maxp;
3223: }
3224: } else {
3225: // Load remote
3226: PetscInt s = 0;
3227: PetscInt *sizes;
3228: MPI_Status status;
3231: MPI_Recv(&min, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
3232: MPI_Recv(&max, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
3233: sieve.setChart(typename ISieve::chart_type(min, max));
3234: PetscMalloc((max - min)*2 * sizeof(PetscInt), &sizes);CHKERRXX(ierr);
3235: MPI_Recv(sizes, (max - min)*2, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
3236: for(typename ISieve::point_type p = min; p < max; ++p, ++s) {
3237: sieve.setConeSize(p, sizes[s*2+0]);
3238: sieve.setSupportSize(p, sizes[s*2+1]);
3239: }
3240: PetscFree(sizes);CHKERRXX(ierr);
3241: }
3242: sieve.allocate();
3243: // Load covering
3244: if (sieve.commRank() == 0) {
3245: // Load local
3246: typename ISieve::index_type maxSize = std::max(sieve.getMaxConeSize(), sieve.getMaxSupportSize());
3247: typename ISieve::point_type *points = new typename ISieve::point_type[maxSize];
3249: for(typename ISieve::point_type p = min; p < max; ++p) {
3250: typename ISieve::index_type coneSize = sieve.getConeSize(p);
3251: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
3253: if (coneSize > 0) {
3254: for(int c = 0; c < coneSize; ++c) {
3255: fs >> points[c];
3256: }
3257: sieve.setCone(points, p);
3258: if (sieve.orientedCones()) {
3259: for(int c = 0; c < coneSize; ++c) {
3260: fs >> points[c];
3261: }
3262: sieve.setConeOrientation(points, p);
3263: }
3264: }
3265: if (supportSize > 0) {
3266: for(int s = 0; s < supportSize; ++s) {
3267: fs >> points[s];
3268: }
3269: sieve.setSupport(p, points);
3270: }
3271: }
3272: delete [] points;
3273: // Load and send remote
3274: for(int pr = 1; pr < sieve.commSize(); ++pr) {
3275: PetscInt size = (maxs[pr] - mins[pr])+totalConeSizes[pr]*2 + (maxs[pr] - mins[pr])+totalSupportSizes[pr];
3276: PetscInt off = 0;
3277: PetscInt *data;
3280: MPI_Send(&size, 1, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
3281: // There is no nice way to make a generic MPI type here. Really sucky
3282: PetscMalloc(size * sizeof(PetscInt), &data);CHKERRXX(ierr);
3283: for(typename ISieve::point_type p = mins[pr]; p < maxs[pr]; ++p) {
3284: PetscInt coneSize, supportSize;
3286: fs >> coneSize;
3287: data[off++] = coneSize;
3288: if (coneSize > 0) {
3289: for(int c = 0; c < coneSize; ++c) {
3290: fs >> data[off++];
3291: }
3292: for(int c = 0; c < coneSize; ++c) {
3293: fs >> data[off++];
3294: }
3295: }
3296: fs >> supportSize;
3297: data[off++] = supportSize;
3298: if (supportSize > 0) {
3299: for(int s = 0; s < supportSize; ++s) {
3300: fs >> data[off++];
3301: }
3302: }
3303: }
3304: assert(off == size);
3305: MPI_Send(data, size, MPIU_INT, pr, 1, sieve.comm());CHKERRXX(ierr);
3306: PetscFree(data);CHKERRXX(ierr);
3307: }
3308: delete [] mins;
3309: delete [] maxs;
3310: } else {
3311: // Load remote
3312: PetscInt size;
3313: PetscInt *data;
3314: PetscInt off = 0;
3315: MPI_Status status;
3318: MPI_Recv(&size, 1, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
3319: PetscMalloc(size*sizeof(PetscInt), &data);CHKERRXX(ierr);
3320: MPI_Recv(data, size, MPIU_INT, 0, 1, sieve.comm(), &status);CHKERRXX(ierr);
3321: for(typename ISieve::point_type p = min; p < max; ++p) {
3322: typename ISieve::index_type coneSize = sieve.getConeSize(p);
3323: typename ISieve::index_type supportSize = sieve.getSupportSize(p);
3324: PetscInt cs = data[off++];
3326: assert(cs == coneSize);
3327: if (coneSize > 0) {
3328: sieve.setCone(&data[off], p);
3329: off += coneSize;
3330: if (sieve.orientedCones()) {
3331: sieve.setConeOrientation(&data[off], p);
3332: off += coneSize;
3333: }
3334: }
3335: PetscInt ss = data[off++];
3337: assert(ss == supportSize);
3338: if (supportSize > 0) {
3339: sieve.setSupport(p, &data[off]);
3340: off += supportSize;
3341: }
3342: }
3343: assert(off == size);
3344: }
3345: // Load renumbering
3346: }
3347: };
3348: }
3350: #endif