Actual source code: Sections.hh
petsc-3.3-p7 2013-05-11
1: #ifndef included_ALE_Sections_hh
2: #define included_ALE_Sections_hh
4: namespace ALE {
5: template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
6: class BaseSection : public ALE::ParallelObject {
7: public:
8: typedef Sieve_ sieve_type;
9: typedef Alloc_ alloc_type;
10: typedef int value_type;
11: typedef typename sieve_type::target_type point_type;
12: typedef typename sieve_type::traits::baseSequence chart_type;
13: protected:
14: Obj<sieve_type> _sieve;
15: chart_type _chart;
16: int _sizes[2];
17: public:
18: BaseSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
19: ~BaseSection() {};
20: public: // Verifiers
21: bool hasPoint(const point_type& point) const {
22: return this->_sieve->baseContains(point);
23: };
24: public:
25: const chart_type& getChart() const {
26: return this->_chart;
27: };
28: int getFiberDimension(const point_type& p) const {
29: return this->hasPoint(p) ? 1 : 0;
30: };
31: const value_type *restrictSpace() const {
32: return this->_sizes;
33: };
34: const value_type *restrictPoint(const point_type& p) const {
35: if (this->hasPoint(p)) return this->_sizes;
36: return &this->_sizes[1];
37: };
38: };
40: template<typename Sieve_, typename Alloc_ = malloc_allocator<int> >
41: class ConeSizeSection : public ALE::ParallelObject {
42: public:
43: typedef Sieve_ sieve_type;
44: typedef Alloc_ alloc_type;
45: typedef int value_type;
46: typedef typename sieve_type::target_type point_type;
47: typedef BaseSection<sieve_type, alloc_type> atlas_type;
48: typedef typename atlas_type::chart_type chart_type;
49: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
50: typedef typename atlas_alloc_type::pointer atlas_ptr;
51: protected:
52: Obj<sieve_type> _sieve;
53: Obj<atlas_type> _atlas;
54: int _size;
55: public:
56: ConeSizeSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
57: atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
58: atlas_alloc_type().construct(pAtlas, atlas_type(sieve));
59: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
60: };
61: ~ConeSizeSection() {};
62: public: // Verifiers
63: bool hasPoint(const point_type& point) {
64: return this->_atlas->hasPoint(point);
65: };
66: public: // Accessors
67: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
68: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
69: public:
70: int getFiberDimension(const point_type& p) {
71: return this->hasPoint(p) ? 1 : 0;
72: };
73: const value_type *restrictPoint(const point_type& p) {
74: this->_size = this->_sieve->cone(p)->size();
75: return &this->_size;
76: };
77: };
79: template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::source_type> >
80: class ConeSection : public ALE::ParallelObject {
81: public:
82: typedef Sieve_ sieve_type;
83: typedef Alloc_ alloc_type;
84: typedef typename sieve_type::target_type point_type;
85: typedef typename sieve_type::source_type value_type;
86: typedef ConeSizeSection<sieve_type, alloc_type> atlas_type;
87: typedef typename atlas_type::chart_type chart_type;
88: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
89: typedef typename atlas_alloc_type::pointer atlas_ptr;
90: protected:
91: Obj<sieve_type> _sieve;
92: Obj<atlas_type> _atlas;
93: alloc_type _allocator;
94: public:
95: ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
96: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
97: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
98: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
99: };
100: ~ConeSection() {};
101: protected:
102: value_type *getRawArray(const int size) {
103: static value_type *array = NULL;
104: static int maxSize = 0;
106: if (size > maxSize) {
107: const value_type dummy(0);
109: if (array) {
110: for(int i = 0; i < maxSize; ++i) {this->_allocator.destroy(array+i);}
111: this->_allocator.deallocate(array, maxSize);
112: }
113: maxSize = size;
114: array = this->_allocator.allocate(maxSize);
115: for(int i = 0; i < maxSize; ++i) {this->_allocator.construct(array+i, dummy);}
116: };
117: return array;
118: };
119: public: // Verifiers
120: bool hasPoint(const point_type& point) {
121: return this->_atlas->hasPoint(point);
122: };
123: public: // Accessors
124: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
125: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
126: public: // Sizes and storage
127: int getFiberDimension(const point_type& p) {
128: return this->_atlas->restrictPoint(p)[0];
129: };
130: public: // Restriction and update
131: const value_type *restrictPoint(const point_type& p) {
132: const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
133: value_type *array = this->getRawArray(cone->size());
134: int c = 0;
136: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
137: array[c++] = *c_iter;
138: }
139: return array;
140: };
141: };
143: template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
144: class BaseSectionV : public ALE::ParallelObject {
145: public:
146: typedef Sieve_ sieve_type;
147: typedef Alloc_ alloc_type;
148: typedef int value_type;
149: typedef typename sieve_type::target_type point_type;
150: //typedef typename sieve_type::traits::baseSequence chart_type;
151: typedef int chart_type;
152: protected:
153: Obj<sieve_type> _sieve;
154: //chart_type _chart;
155: int _sizes[2];
156: public:
157: //BaseSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
158: BaseSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {_sizes[0] = 1; _sizes[1] = 0;};
159: ~BaseSectionV() {};
160: public: // Verifiers
161: bool hasPoint(const point_type& point) const {
162: return this->_sieve->baseContains(point);
163: };
164: public:
165: //const chart_type& getChart() const {
166: // return this->_chart;
167: //};
168: int getFiberDimension(const point_type& p) const {
169: return this->hasPoint(p) ? 1 : 0;
170: };
171: const value_type *restrictSpace() const {
172: return this->_sizes;
173: };
174: const value_type *restrictPoint(const point_type& p) const {
175: if (this->hasPoint(p)) return this->_sizes;
176: return &this->_sizes[1];
177: };
178: };
180: template<typename Sieve_, typename Alloc_ = malloc_allocator<int> >
181: class ConeSizeSectionV : public ALE::ParallelObject {
182: public:
183: typedef Sieve_ sieve_type;
184: typedef Alloc_ alloc_type;
185: typedef int value_type;
186: typedef typename sieve_type::target_type point_type;
187: typedef BaseSectionV<sieve_type, alloc_type> atlas_type;
188: typedef typename atlas_type::chart_type chart_type;
189: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
190: typedef typename atlas_alloc_type::pointer atlas_ptr;
191: protected:
192: Obj<sieve_type> _sieve;
193: Obj<atlas_type> _atlas;
194: int _size;
195: public:
196: ConeSizeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
197: atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
198: atlas_alloc_type().construct(pAtlas, atlas_type(sieve));
199: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
200: };
201: ~ConeSizeSectionV() {};
202: public: // Verifiers
203: bool hasPoint(const point_type& point) {
204: return this->_atlas->hasPoint(point);
205: };
206: public: // Accessors
207: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
208: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
209: public:
210: int getFiberDimension(const point_type& p) {
211: return this->hasPoint(p) ? 1 : 0;
212: };
213: const value_type *restrictPoint(const point_type& p) {
214: this->_size = this->_sieve->getConeSize(p);
215: return &this->_size;
216: };
217: };
219: template<typename Sieve_, typename Alloc_ = malloc_allocator<typename Sieve_::source_type> >
220: class ConeSectionV : public ALE::ParallelObject {
221: public:
222: typedef Sieve_ sieve_type;
223: typedef Alloc_ alloc_type;
224: typedef typename sieve_type::target_type point_type;
225: typedef typename sieve_type::source_type value_type;
226: typedef ConeSizeSectionV<sieve_type, alloc_type> atlas_type;
227: typedef typename atlas_type::chart_type chart_type;
228: typedef typename ISieveVisitor::PointRetriever<sieve_type> visitor_type;
229: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
230: typedef typename atlas_alloc_type::pointer atlas_ptr;
231: protected:
232: Obj<sieve_type> _sieve;
233: Obj<atlas_type> _atlas;
234: visitor_type *_cV;
235: alloc_type _allocator;
236: public:
237: ConeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
238: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
239: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
240: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
241: this->_cV = new visitor_type(std::max(0, sieve->getMaxConeSize()));
242: };
243: ~ConeSectionV() {
244: delete this->_cV;
245: };
246: public: // Verifiers
247: bool hasPoint(const point_type& point) {
248: return this->_atlas->hasPoint(point);
249: };
250: public: // Accessors
251: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
252: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
253: public: // Sizes and storage
254: int getFiberDimension(const point_type& p) {
255: return this->_atlas->restrictPoint(p)[0];
256: };
257: public: // Restriction and update
258: const value_type *restrictPoint(const point_type& p) {
259: this->_cV->clear();
260: this->_sieve->cone(p, *this->_cV);
261: return this->_cV->getPoints();
262: };
263: };
265: template<typename Sieve_, typename Alloc_ = malloc_allocator<OrientedPoint<typename Sieve_::source_type> > >
266: class OrientedConeSectionV : public ALE::ParallelObject {
267: public:
268: typedef Sieve_ sieve_type;
269: typedef Alloc_ alloc_type;
270: typedef typename sieve_type::target_type point_type;
271: typedef OrientedPoint<typename sieve_type::source_type> value_type;
272: typedef typename alloc_type::template rebind<int>::other int_alloc_type;
273: typedef ConeSizeSectionV<sieve_type, int_alloc_type> atlas_type;
274: typedef typename atlas_type::chart_type chart_type;
275: typedef typename ISieveVisitor::PointRetriever<sieve_type> visitor_type;
276: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
277: typedef typename atlas_alloc_type::pointer atlas_ptr;
278: protected:
279: Obj<sieve_type> _sieve;
280: Obj<atlas_type> _atlas;
281: visitor_type *_cV;
282: alloc_type _allocator;
283: public:
284: OrientedConeSectionV(const Obj<sieve_type>& sieve) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve) {
285: atlas_ptr pAtlas = atlas_alloc_type(this->_allocator).allocate(1);
286: atlas_alloc_type(this->_allocator).construct(pAtlas, atlas_type(sieve));
287: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
288: this->_cV = new visitor_type(std::max(0, sieve->getMaxConeSize()));
289: };
290: ~OrientedConeSectionV() {
291: delete this->_cV;
292: };
293: public: // Verifiers
294: bool hasPoint(const point_type& point) {
295: return this->_atlas->hasPoint(point);
296: };
297: public: // Accessors
298: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
299: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
300: public: // Sizes and storage
301: int getFiberDimension(const point_type& p) {
302: return this->_atlas->restrictPoint(p)[0];
303: };
304: int size() {
305: int s = 0;
306: for(int p = this->_sieve->getChart().min(); p < this->_sieve->getChart().max(); ++p) {
307: s += this->getFiberDimension(p);
308: }
309: return s;
310: };
311: public: // Restriction and update
312: const value_type *restrictPoint(const point_type& p) {
313: this->_cV->clear();
314: this->_sieve->orientedCone(p, *this->_cV);
315: return (const value_type *) this->_cV->getOrientedPoints();
316: };
317: };
319: template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
320: class LabelBaseSection : public ALE::ParallelObject {
321: public:
322: typedef Sieve_ sieve_type;
323: typedef Label_ label_type;
324: typedef Alloc_ alloc_type;
325: typedef int value_type;
326: typedef typename sieve_type::target_type point_type;
327: typedef typename sieve_type::traits::baseSequence chart_type;
328: protected:
329: Obj<sieve_type> _sieve;
330: Obj<label_type> _label;
331: chart_type _chart;
332: int _sizes[2];
333: public:
334: LabelBaseSection(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
335: ~LabelBaseSection() {};
336: public: // Verifiers
337: bool hasPoint(const point_type& point) const {
338: return this->_label->cone(point)->size() ? true : false;
339: };
340: public:
341: const chart_type& getChart() const {
342: return this->_chart;
343: };
344: int getFiberDimension(const point_type& p) const {
345: return this->hasPoint(p) ? 1 : 0;
346: };
347: const value_type *restrictSpace() const {
348: return this->_sizes;
349: };
350: const value_type *restrictPoint(const point_type& p) const {
351: if (this->hasPoint(p)) return this->_sizes;
352: return &this->_sizes[1];
353: };
354: };
356: template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<int>, typename Atlas_ = LabelBaseSection<Sieve_, Label_, Alloc_> >
357: class LabelSection : public ALE::ParallelObject {
358: public:
359: typedef Sieve_ sieve_type;
360: typedef Label_ label_type;
361: typedef Alloc_ alloc_type;
362: typedef int value_type;
363: typedef typename sieve_type::target_type point_type;
364: typedef Atlas_ atlas_type;
365: typedef typename atlas_type::chart_type chart_type;
366: typedef typename alloc_type::template rebind<atlas_type>::other atlas_alloc_type;
367: typedef typename atlas_alloc_type::pointer atlas_ptr;
368: protected:
369: Obj<sieve_type> _sieve;
370: Obj<label_type> _label;
371: Obj<atlas_type> _atlas;
372: int _size;
373: int _value;
374: public:
375: LabelSection(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label) {
376: atlas_ptr pAtlas = atlas_alloc_type().allocate(1);
377: atlas_alloc_type().construct(pAtlas, atlas_type(sieve, label));
378: this->_atlas = Obj<atlas_type>(pAtlas, sizeof(atlas_type));
379: };
380: ~LabelSection() {};
381: public: // Verifiers
382: bool hasPoint(const point_type& point) {
383: return this->_atlas->hasPoint(point);
384: };
385: public: // Accessors
386: const Obj<atlas_type>& getAtlas() {return this->_atlas;};
387: void setAtlas(const Obj<atlas_type>& atlas) {this->_atlas = atlas;};
388: public:
389: int getFiberDimension(const point_type& p) {
390: return this->hasPoint(p) ? 1 : 0;
391: };
392: const value_type *restrictPoint(const point_type& p) {
393: this->_value = *this->_label->cone(p)->begin();
394: return &this->_value;
395: };
396: };
398: template<typename Sieve_, typename Label_, typename Alloc_ = malloc_allocator<typename Sieve_::target_type> >
399: class LabelBaseSectionV : public ALE::ParallelObject {
400: public:
401: typedef Sieve_ sieve_type;
402: typedef Label_ label_type;
403: typedef Alloc_ alloc_type;
404: typedef int value_type;
405: typedef typename sieve_type::target_type point_type;
406: //typedef typename sieve_type::traits::baseSequence chart_type;
407: typedef int chart_type;
408: protected:
409: Obj<sieve_type> _sieve;
410: Obj<label_type> _label;
411: //chart_type _chart;
412: int _sizes[2];
413: public:
414: //LabelBaseSectionV(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label), _chart(*sieve->base()) {_sizes[0] = 1; _sizes[1] = 0;};
415: LabelBaseSectionV(const Obj<sieve_type>& sieve, const Obj<label_type>& label) : ParallelObject(sieve->comm(), sieve->debug()), _sieve(sieve), _label(label) {_sizes[0] = 1; _sizes[1] = 0;};
416: ~LabelBaseSectionV() {};
417: public: // Verifiers
418: bool hasPoint(const point_type& point) const {
419: return this->_label->cone(point)->size() ? true : false;
420: };
421: public:
422: //const chart_type& getChart() const {
423: // return this->_chart;
424: //};
425: int getFiberDimension(const point_type& p) const {
426: return this->hasPoint(p) ? 1 : 0;
427: };
428: const value_type *restrictSpace() const {
429: return this->_sizes;
430: };
431: const value_type *restrictPoint(const point_type& p) const {
432: if (this->hasPoint(p)) return this->_sizes;
433: return &this->_sizes[1];
434: };
435: };
437: namespace New {
438: // This section takes an existing section, and reports instead the fiber dimensions as values
439: template<typename Section_>
440: class SizeSection : public ALE::ParallelObject {
441: public:
442: typedef Section_ section_type;
443: typedef typename section_type::point_type point_type;
444: typedef int value_type;
445: protected:
446: Obj<section_type> _section;
447: value_type _size;
448: public:
449: SizeSection(const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, section->debug()), _section(section) {};
450: virtual ~SizeSection() {};
451: public:
452: bool hasPoint(const point_type& point) {
453: return this->_section->hasPoint(point);
454: };
455: const value_type *restrictPoint(const point_type& p) {
456: this->_size = this->_section->getFiberDimension(p);
457: return &this->_size;
458: };
459: public:
460: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
461: this->_section->view(name, comm);
462: };
463: };
465: // This section reports as values the size of the partition associated with the partition point
466: template<typename Bundle_, typename Marker_>
467: class PartitionSizeSection : public ALE::ParallelObject {
468: public:
469: typedef Bundle_ bundle_type;
470: typedef typename bundle_type::sieve_type sieve_type;
471: typedef typename bundle_type::point_type point_type;
472: typedef Marker_ marker_type;
473: typedef int value_type;
474: typedef std::map<marker_type, int> sizes_type;
475: protected:
476: sizes_type _sizes;
477: int _height;
478: void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
479: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(this->_height);
480: const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
481: std::map<marker_type, std::set<point_type> > points;
483: if (numElements != (int) cells->size()) {
484: throw ALE::Exception("Partition size does not match the number of elements");
485: }
486: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
487: typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
488: const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
489: const int idx = cNumbering->getIndex(*e_iter);
491: points[partition[idx]].insert(closure->begin(), closure->end());
492: if (this->_height > 0) {
493: const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);
495: points[partition[idx]].insert(star->begin(), star->end());
496: }
497: }
498: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
499: this->_sizes[p_iter->first] = p_iter->second.size();
500: }
501: };
502: public:
503: PartitionSizeSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
504: this->_init(bundle, numElements, partition);
505: };
506: virtual ~PartitionSizeSection() {};
507: public:
508: bool hasPoint(const point_type& point) {return true;};
509: const value_type *restrictPoint(const point_type& p) {
510: return &this->_sizes[p];
511: };
512: public:
513: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
514: ostringstream txt;
515: int rank;
517: if (comm == MPI_COMM_NULL) {
518: comm = this->comm();
519: rank = this->commRank();
520: } else {
521: MPI_Comm_rank(comm, &rank);
522: }
523: if (name == "") {
524: if(rank == 0) {
525: txt << "viewing a PartitionSizeSection" << std::endl;
526: }
527: } else {
528: if(rank == 0) {
529: txt << "viewing PartitionSizeSection '" << name << "'" << std::endl;
530: }
531: }
532: for(typename sizes_type::const_iterator s_iter = this->_sizes.begin(); s_iter != this->_sizes.end(); ++s_iter) {
533: const marker_type& partition = s_iter->first;
534: const value_type size = s_iter->second;
536: txt << "[" << this->commRank() << "]: Partition " << partition << " size " << size << std::endl;
537: }
538: PetscSynchronizedPrintf(comm, txt.str().c_str());
539: PetscSynchronizedFlush(comm);
540: };
541: };
543: template<typename Point_>
544: class PartitionDomain {
545: public:
546: typedef Point_ point_type;
547: public:
548: PartitionDomain() {};
549: ~PartitionDomain() {};
550: public:
551: int count(const point_type& point) const {return 1;};
552: };
554: // This section returns the points in each partition
555: template<typename Bundle_, typename Marker_>
556: class PartitionSection : public ALE::ParallelObject {
557: public:
558: typedef Bundle_ bundle_type;
559: typedef typename bundle_type::sieve_type sieve_type;
560: typedef typename bundle_type::point_type point_type;
561: typedef Marker_ marker_type;
562: typedef int value_type;
563: typedef std::map<marker_type, point_type*> points_type;
564: typedef PartitionDomain<point_type> chart_type;
565: protected:
566: points_type _points;
567: chart_type _domain;
568: int _height;
569: void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
570: // Should check for patch 0
571: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(this->_height);
572: const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
573: std::map<marker_type, std::set<point_type> > points;
574: std::map<marker_type, int> offsets;
576: if (numElements != (int) cells->size()) {
577: throw ALE::Exception("Partition size does not match the number of elements");
578: }
579: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
580: typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
581: const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
582: const int idx = cNumbering->getIndex(*e_iter);
584: points[partition[idx]].insert(closure->begin(), closure->end());
585: if (this->_height > 0) {
586: const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);
588: points[partition[idx]].insert(star->begin(), star->end());
589: }
590: }
591: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
592: this->_points[p_iter->first] = new point_type[p_iter->second.size()];
593: offsets[p_iter->first] = 0;
594: for(typename std::set<point_type>::const_iterator s_iter = p_iter->second.begin(); s_iter != p_iter->second.end(); ++s_iter) {
595: this->_points[p_iter->first][offsets[p_iter->first]++] = *s_iter;
596: }
597: }
598: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
599: if (offsets[p_iter->first] != (int) p_iter->second.size()) {
600: ostringstream txt;
601: txt << "Invalid offset for partition " << p_iter->first << ": " << offsets[p_iter->first] << " should be " << p_iter->second.size();
602: throw ALE::Exception(txt.str().c_str());
603: }
604: }
605: };
606: public:
607: PartitionSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
608: this->_init(bundle, numElements, partition);
609: };
610: virtual ~PartitionSection() {
611: for(typename points_type::iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
612: delete [] p_iter->second;
613: }
614: };
615: public:
616: const chart_type& getChart() {return this->_domain;};
617: bool hasPoint(const point_type& point) {return true;};
618: const value_type *restrictPoint(const point_type& p) {
619: return this->_points[p];
620: };
621: public:
622: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
623: ostringstream txt;
624: int rank;
626: if (comm == MPI_COMM_NULL) {
627: comm = this->comm();
628: rank = this->commRank();
629: } else {
630: MPI_Comm_rank(comm, &rank);
631: }
632: if (name == "") {
633: if(rank == 0) {
634: txt << "viewing a PartitionSection" << std::endl;
635: }
636: } else {
637: if(rank == 0) {
638: txt << "viewing PartitionSection '" << name << "'" << std::endl;
639: }
640: }
641: for(typename points_type::const_iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
642: const marker_type& partition = p_iter->first;
643: //const point_type *points = p_iter->second;
645: txt << "[" << this->commRank() << "]: Partition " << partition << std::endl;
646: }
647: if (this->_points.size() == 0) {
648: txt << "[" << this->commRank() << "]: empty" << std::endl;
649: }
650: PetscSynchronizedPrintf(comm, txt.str().c_str());
651: PetscSynchronizedFlush(comm);
652: };
653: };
655: template<typename Bundle_, typename Sieve_>
656: class ConeSizeSection : public ALE::ParallelObject {
657: public:
658: typedef ConeSizeSection<Bundle_, Sieve_> section_type;
659: typedef int patch_type;
660: typedef Bundle_ bundle_type;
661: typedef Sieve_ sieve_type;
662: typedef typename bundle_type::point_type point_type;
663: typedef int value_type;
664: protected:
665: Obj<bundle_type> _bundle;
666: Obj<sieve_type> _sieve;
667: value_type _size;
668: int _minHeight;
669: Obj<section_type> _section;
670: public:
671: ConeSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumHeight = 0) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _bundle(bundle), _sieve(sieve), _minHeight(minimumHeight) {
672: this->_section = this;
673: this->_section.addRef();
674: };
675: virtual ~ConeSizeSection() {};
676: public: // Verifiers
677: bool hasPoint(const point_type& point) {return true;};
678: public: // Restriction
679: const value_type *restrictPoint(const point_type& p) {
680: if ((this->_minHeight == 0) || (this->_bundle->height(p) >= this->_minHeight)) {
681: this->_size = this->_sieve->cone(p)->size();
682: } else {
683: this->_size = 0;
684: }
685: return &this->_size;
686: };
687: public: // Adapter
688: const Obj<section_type>& getSection(const patch_type& patch) {
689: return this->_section;
690: };
691: public:
692: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
693: ostringstream txt;
694: int rank;
696: if (comm == MPI_COMM_NULL) {
697: comm = this->comm();
698: rank = this->commRank();
699: } else {
700: MPI_Comm_rank(comm, &rank);
701: }
702: if (name == "") {
703: if(rank == 0) {
704: txt << "viewing a ConeSizeSection" << std::endl;
705: }
706: } else {
707: if(rank == 0) {
708: txt << "viewing ConeSizeSection '" << name << "'" << std::endl;
709: }
710: }
711: PetscSynchronizedPrintf(comm, txt.str().c_str());
712: PetscSynchronizedFlush(comm);
713: };
714: };
716: template<typename Sieve_>
717: class ConeSection : public ALE::ParallelObject {
718: public:
719: typedef Sieve_ sieve_type;
720: typedef typename sieve_type::target_type point_type;
721: typedef typename sieve_type::source_type value_type;
722: typedef PartitionDomain<sieve_type> chart_type;
723: protected:
724: Obj<sieve_type> _sieve;
725: int _coneSize;
726: value_type *_cone;
727: chart_type _domain;
728: void ensureCone(const int size) {
729: if (size > this->_coneSize) {
730: if (this->_cone) delete [] this->_cone;
731: this->_coneSize = size;
732: this->_cone = new value_type[this->_coneSize];
733: }
734: };
735: public:
736: ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _coneSize(-1), _cone(NULL) {};
737: virtual ~ConeSection() {if (this->_cone) delete [] this->_cone;};
738: public:
739: const chart_type& getChart() {return this->_domain;};
740: bool hasPoint(const point_type& point) {return true;};
741: const value_type *restrictPoint(const point_type& p) {
742: const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
743: int c = 0;
745: this->ensureCone(cone->size());
746: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
747: this->_cone[c++] = *c_iter;
748: }
749: return this->_cone;
750: };
751: public:
752: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
753: ostringstream txt;
754: int rank;
756: if (comm == MPI_COMM_NULL) {
757: comm = this->comm();
758: rank = this->commRank();
759: } else {
760: MPI_Comm_rank(comm, &rank);
761: }
762: if (name == "") {
763: if(rank == 0) {
764: txt << "viewing a ConeSection" << std::endl;
765: }
766: } else {
767: if(rank == 0) {
768: txt << "viewing ConeSection '" << name << "'" << std::endl;
769: }
770: }
771: PetscSynchronizedPrintf(comm, txt.str().c_str());
772: PetscSynchronizedFlush(comm);
773: };
774: };
776: template<typename Bundle_, typename Sieve_>
777: class SupportSizeSection : public ALE::ParallelObject {
778: public:
779: typedef Bundle_ bundle_type;
780: typedef Sieve_ sieve_type;
781: typedef typename sieve_type::source_type point_type;
782: typedef typename sieve_type::target_type value_type;
783: protected:
784: Obj<bundle_type> _bundle;
785: Obj<sieve_type> _sieve;
786: value_type _size;
787: int _minDepth;
788: public:
789: SupportSizeSection(const Obj<bundle_type>& bundle, const Obj<sieve_type>& sieve, int minimumDepth = 0) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _bundle(bundle), _sieve(sieve), _minDepth(minimumDepth) {};
790: virtual ~SupportSizeSection() {};
791: public:
792: bool hasPoint(const point_type& point) {return true;};
793: const value_type *restrictPoint(const point_type& p) {
794: if ((this->_minDepth == 0) || (this->_bundle->depth(p) >= this->_minDepth)) {
795: this->_size = this->_sieve->support(p)->size();
796: } else {
797: this->_size = 0;
798: }
799: return &this->_size;
800: };
801: public:
802: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
803: ostringstream txt;
804: int rank;
806: if (comm == MPI_COMM_NULL) {
807: comm = this->comm();
808: rank = this->commRank();
809: } else {
810: MPI_Comm_rank(comm, &rank);
811: }
812: if (name == "") {
813: if(rank == 0) {
814: txt << "viewing a SupportSizeSection" << std::endl;
815: }
816: } else {
817: if(rank == 0) {
818: txt << "viewing SupportSizeSection '" << name << "'" << std::endl;
819: }
820: }
821: PetscSynchronizedPrintf(comm, txt.str().c_str());
822: PetscSynchronizedFlush(comm);
823: };
824: };
826: template<typename Sieve_>
827: class SupportSection : public ALE::ParallelObject {
828: public:
829: typedef Sieve_ sieve_type;
830: typedef typename sieve_type::source_type point_type;
831: typedef typename sieve_type::target_type value_type;
832: typedef PartitionDomain<sieve_type> chart_type;
833: protected:
834: Obj<sieve_type> _sieve;
835: int _supportSize;
836: value_type *_support;
837: chart_type _domain;
838: void ensureSupport(const int size) {
839: if (size > this->_supportSize) {
840: if (this->_support) delete [] this->_support;
841: this->_supportSize = size;
842: this->_support = new value_type[this->_supportSize];
843: }
844: };
845: public:
846: SupportSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _supportSize(-1), _support(NULL) {};
847: virtual ~SupportSection() {if (this->_support) delete [] this->_support;};
848: public:
849: const chart_type& getChart() {return this->_domain;};
850: bool hasPoint(const point_type& point) {return true;};
851: const value_type *restrictPoint(const point_type& p) {
852: const Obj<typename sieve_type::traits::supportSequence>& support = this->_sieve->support(p);
853: int s = 0;
855: this->ensureSupport(support->size());
856: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
857: this->_support[s++] = *s_iter;
858: }
859: return this->_support;
860: };
861: public:
862: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
863: ostringstream txt;
864: int rank;
866: if (comm == MPI_COMM_NULL) {
867: comm = this->comm();
868: rank = this->commRank();
869: } else {
870: MPI_Comm_rank(comm, &rank);
871: }
872: if (name == "") {
873: if(rank == 0) {
874: txt << "viewing a SupportSection" << std::endl;
875: }
876: } else {
877: if(rank == 0) {
878: txt << "viewing SupportSection '" << name << "'" << std::endl;
879: }
880: }
881: PetscSynchronizedPrintf(comm, txt.str().c_str());
882: PetscSynchronizedFlush(comm);
883: };
884: };
886: template<typename Sieve_, typename Section_>
887: class ArrowSection : public ALE::ParallelObject {
888: public:
889: typedef Sieve_ sieve_type;
890: typedef Section_ section_type;
891: typedef typename sieve_type::target_type point_type;
892: typedef typename section_type::point_type arrow_type;
893: typedef typename section_type::value_type value_type;
894: protected:
895: Obj<sieve_type> _sieve;
896: Obj<section_type> _section;
897: int _coneSize;
898: value_type *_cone;
899: void ensureCone(const int size) {
900: if (size > this->_coneSize) {
901: if (this->_cone) delete [] this->_cone;
902: this->_coneSize = size;
903: this->_cone = new value_type[this->_coneSize];
904: }
905: };
906: public:
907: ArrowSection(const Obj<sieve_type>& sieve, const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _section(section), _coneSize(-1), _cone(NULL) {};
908: virtual ~ArrowSection() {if (this->_cone) delete [] this->_cone;};
909: public:
910: bool hasPoint(const point_type& point) {return this->_sieve->baseContains(point);};
911: const value_type *restrictPoint(const point_type& p) {
912: const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
913: int c = -1;
915: this->ensureCone(cone->size());
916: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
917: this->_cone[++c] = this->_section->restrictPoint(arrow_type(*c_iter, p))[0];
918: }
919: return this->_cone;
920: };
921: public:
922: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
923: ostringstream txt;
924: int rank;
926: if (comm == MPI_COMM_NULL) {
927: comm = this->comm();
928: rank = this->commRank();
929: } else {
930: MPI_Comm_rank(comm, &rank);
931: }
932: if (name == "") {
933: if(rank == 0) {
934: txt << "viewing a ConeSection" << std::endl;
935: }
936: } else {
937: if(rank == 0) {
938: txt << "viewing ConeSection '" << name << "'" << std::endl;
939: }
940: }
941: PetscSynchronizedPrintf(comm, txt.str().c_str());
942: PetscSynchronizedFlush(comm);
943: };
944: };
945: }
946: }
947: #endif