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