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