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> >
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 LabelBaseSection<sieve_type, label_type, alloc_type> 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: namespace New {
396: // This section takes an existing section, and reports instead the fiber dimensions as values
397: template<typename Section_>
398: class SizeSection : public ALE::ParallelObject {
399: public:
400: typedef Section_ section_type;
401: typedef typename section_type::point_type point_type;
402: typedef int value_type;
403: protected:
404: Obj<section_type> _section;
405: value_type _size;
406: public:
407: SizeSection(const Obj<section_type>& section) : ParallelObject(MPI_COMM_SELF, section->debug()), _section(section) {};
408: virtual ~SizeSection() {};
409: public:
410: bool hasPoint(const point_type& point) {
411: return this->_section->hasPoint(point);
412: };
413: const value_type *restrictPoint(const point_type& p) {
414: this->_size = this->_section->getFiberDimension(p);
415: return &this->_size;
416: };
417: public:
418: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
419: this->_section->view(name, comm);
420: };
421: };
423: // This section reports as values the size of the partition associated with the partition point
424: template<typename Bundle_, typename Marker_>
425: class PartitionSizeSection : public ALE::ParallelObject {
426: public:
427: typedef Bundle_ bundle_type;
428: typedef typename bundle_type::sieve_type sieve_type;
429: typedef typename bundle_type::point_type point_type;
430: typedef Marker_ marker_type;
431: typedef int value_type;
432: typedef std::map<marker_type, int> sizes_type;
433: protected:
434: sizes_type _sizes;
435: int _height;
436: void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
437: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(this->_height);
438: const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
439: std::map<marker_type, std::set<point_type> > points;
441: if (numElements != (int) cells->size()) {
442: throw ALE::Exception("Partition size does not match the number of elements");
443: }
444: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
445: typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
446: const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
447: const int idx = cNumbering->getIndex(*e_iter);
449: points[partition[idx]].insert(closure->begin(), closure->end());
450: if (this->_height > 0) {
451: const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);
453: points[partition[idx]].insert(star->begin(), star->end());
454: }
455: }
456: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
457: this->_sizes[p_iter->first] = p_iter->second.size();
458: }
459: };
460: public:
461: PartitionSizeSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
462: this->_init(bundle, numElements, partition);
463: };
464: virtual ~PartitionSizeSection() {};
465: public:
466: bool hasPoint(const point_type& point) {return true;};
467: const value_type *restrictPoint(const point_type& p) {
468: return &this->_sizes[p];
469: };
470: public:
471: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
472: ostringstream txt;
473: int rank;
475: if (comm == MPI_COMM_NULL) {
476: comm = this->comm();
477: rank = this->commRank();
478: } else {
479: MPI_Comm_rank(comm, &rank);
480: }
481: if (name == "") {
482: if(rank == 0) {
483: txt << "viewing a PartitionSizeSection" << std::endl;
484: }
485: } else {
486: if(rank == 0) {
487: txt << "viewing PartitionSizeSection '" << name << "'" << std::endl;
488: }
489: }
490: for(typename sizes_type::const_iterator s_iter = this->_sizes.begin(); s_iter != this->_sizes.end(); ++s_iter) {
491: const marker_type& partition = s_iter->first;
492: const value_type size = s_iter->second;
494: txt << "[" << this->commRank() << "]: Partition " << partition << " size " << size << std::endl;
495: }
496: PetscSynchronizedPrintf(comm, txt.str().c_str());
497: PetscSynchronizedFlush(comm);
498: };
499: };
501: template<typename Point_>
502: class PartitionDomain {
503: public:
504: typedef Point_ point_type;
505: public:
506: PartitionDomain() {};
507: ~PartitionDomain() {};
508: public:
509: int count(const point_type& point) const {return 1;};
510: };
512: // This section returns the points in each partition
513: template<typename Bundle_, typename Marker_>
514: class PartitionSection : public ALE::ParallelObject {
515: public:
516: typedef Bundle_ bundle_type;
517: typedef typename bundle_type::sieve_type sieve_type;
518: typedef typename bundle_type::point_type point_type;
519: typedef Marker_ marker_type;
520: typedef int value_type;
521: typedef std::map<marker_type, point_type*> points_type;
522: typedef PartitionDomain<point_type> chart_type;
523: protected:
524: points_type _points;
525: chart_type _domain;
526: int _height;
527: void _init(const Obj<bundle_type>& bundle, const int numElements, const marker_type partition[]) {
528: // Should check for patch 0
529: const Obj<typename bundle_type::label_sequence>& cells = bundle->heightStratum(this->_height);
530: const Obj<typename bundle_type::numbering_type>& cNumbering = bundle->getFactory()->getLocalNumbering(bundle, bundle->depth() - this->_height);
531: std::map<marker_type, std::set<point_type> > points;
532: std::map<marker_type, int> offsets;
534: if (numElements != (int) cells->size()) {
535: throw ALE::Exception("Partition size does not match the number of elements");
536: }
537: for(typename bundle_type::label_sequence::iterator e_iter = cells->begin(); e_iter != cells->end(); ++e_iter) {
538: typedef ALE::SieveAlg<bundle_type> sieve_alg_type;
539: const Obj<typename sieve_alg_type::coneArray>& closure = sieve_alg_type::closure(bundle, *e_iter);
540: const int idx = cNumbering->getIndex(*e_iter);
542: points[partition[idx]].insert(closure->begin(), closure->end());
543: if (this->_height > 0) {
544: const Obj<typename sieve_alg_type::supportArray>& star = sieve_alg_type::star(bundle, *e_iter);
546: points[partition[idx]].insert(star->begin(), star->end());
547: }
548: }
549: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
550: this->_points[p_iter->first] = new point_type[p_iter->second.size()];
551: offsets[p_iter->first] = 0;
552: for(typename std::set<point_type>::const_iterator s_iter = p_iter->second.begin(); s_iter != p_iter->second.end(); ++s_iter) {
553: this->_points[p_iter->first][offsets[p_iter->first]++] = *s_iter;
554: }
555: }
556: for(typename std::map<marker_type, std::set<point_type> >::const_iterator p_iter = points.begin(); p_iter != points.end(); ++p_iter) {
557: if (offsets[p_iter->first] != (int) p_iter->second.size()) {
558: ostringstream txt;
559: txt << "Invalid offset for partition " << p_iter->first << ": " << offsets[p_iter->first] << " should be " << p_iter->second.size();
560: throw ALE::Exception(txt.str().c_str());
561: }
562: }
563: };
564: public:
565: PartitionSection(const Obj<bundle_type>& bundle, const int elementHeight, const int numElements, const marker_type *partition) : ParallelObject(MPI_COMM_SELF, bundle->debug()), _height(elementHeight) {
566: this->_init(bundle, numElements, partition);
567: };
568: virtual ~PartitionSection() {
569: for(typename points_type::iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
570: delete [] p_iter->second;
571: }
572: };
573: public:
574: const chart_type& getChart() {return this->_domain;};
575: bool hasPoint(const point_type& point) {return true;};
576: const value_type *restrictPoint(const point_type& p) {
577: return this->_points[p];
578: };
579: public:
580: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
581: ostringstream txt;
582: int rank;
584: if (comm == MPI_COMM_NULL) {
585: comm = this->comm();
586: rank = this->commRank();
587: } else {
588: MPI_Comm_rank(comm, &rank);
589: }
590: if (name == "") {
591: if(rank == 0) {
592: txt << "viewing a PartitionSection" << std::endl;
593: }
594: } else {
595: if(rank == 0) {
596: txt << "viewing PartitionSection '" << name << "'" << std::endl;
597: }
598: }
599: for(typename points_type::const_iterator p_iter = this->_points.begin(); p_iter != this->_points.end(); ++p_iter) {
600: const marker_type& partition = p_iter->first;
601: //const point_type *points = p_iter->second;
603: txt << "[" << this->commRank() << "]: Partition " << partition << std::endl;
604: }
605: if (this->_points.size() == 0) {
606: txt << "[" << this->commRank() << "]: empty" << std::endl;
607: }
608: PetscSynchronizedPrintf(comm, txt.str().c_str());
609: PetscSynchronizedFlush(comm);
610: };
611: };
613: template<typename Bundle_, typename Sieve_>
614: class ConeSizeSection : public ALE::ParallelObject {
615: public:
616: typedef ConeSizeSection<Bundle_, Sieve_> section_type;
617: typedef int patch_type;
618: typedef Bundle_ bundle_type;
619: typedef Sieve_ sieve_type;
620: typedef typename bundle_type::point_type point_type;
621: typedef int value_type;
622: protected:
623: Obj<bundle_type> _bundle;
624: Obj<sieve_type> _sieve;
625: value_type _size;
626: int _minHeight;
627: Obj<section_type> _section;
628: public:
629: 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) {
630: this->_section = this;
631: this->_section.addRef();
632: };
633: virtual ~ConeSizeSection() {};
634: public: // Verifiers
635: bool hasPoint(const point_type& point) {return true;};
636: public: // Restriction
637: const value_type *restrictPoint(const point_type& p) {
638: if ((this->_minHeight == 0) || (this->_bundle->height(p) >= this->_minHeight)) {
639: this->_size = this->_sieve->cone(p)->size();
640: } else {
641: this->_size = 0;
642: }
643: return &this->_size;
644: };
645: public: // Adapter
646: const Obj<section_type>& getSection(const patch_type& patch) {
647: return this->_section;
648: };
649: public:
650: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
651: ostringstream txt;
652: int rank;
654: if (comm == MPI_COMM_NULL) {
655: comm = this->comm();
656: rank = this->commRank();
657: } else {
658: MPI_Comm_rank(comm, &rank);
659: }
660: if (name == "") {
661: if(rank == 0) {
662: txt << "viewing a ConeSizeSection" << std::endl;
663: }
664: } else {
665: if(rank == 0) {
666: txt << "viewing ConeSizeSection '" << name << "'" << std::endl;
667: }
668: }
669: PetscSynchronizedPrintf(comm, txt.str().c_str());
670: PetscSynchronizedFlush(comm);
671: };
672: };
674: template<typename Sieve_>
675: class ConeSection : public ALE::ParallelObject {
676: public:
677: typedef Sieve_ sieve_type;
678: typedef typename sieve_type::target_type point_type;
679: typedef typename sieve_type::source_type value_type;
680: typedef PartitionDomain<sieve_type> chart_type;
681: protected:
682: Obj<sieve_type> _sieve;
683: int _coneSize;
684: value_type *_cone;
685: chart_type _domain;
686: void ensureCone(const int size) {
687: if (size > this->_coneSize) {
688: if (this->_cone) delete [] this->_cone;
689: this->_coneSize = size;
690: this->_cone = new value_type[this->_coneSize];
691: }
692: };
693: public:
694: ConeSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _coneSize(-1), _cone(NULL) {};
695: virtual ~ConeSection() {if (this->_cone) delete [] this->_cone;};
696: public:
697: const chart_type& getChart() {return this->_domain;};
698: bool hasPoint(const point_type& point) {return true;};
699: const value_type *restrictPoint(const point_type& p) {
700: const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
701: int c = 0;
703: this->ensureCone(cone->size());
704: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
705: this->_cone[c++] = *c_iter;
706: }
707: return this->_cone;
708: };
709: public:
710: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
711: ostringstream txt;
712: int rank;
714: if (comm == MPI_COMM_NULL) {
715: comm = this->comm();
716: rank = this->commRank();
717: } else {
718: MPI_Comm_rank(comm, &rank);
719: }
720: if (name == "") {
721: if(rank == 0) {
722: txt << "viewing a ConeSection" << std::endl;
723: }
724: } else {
725: if(rank == 0) {
726: txt << "viewing ConeSection '" << name << "'" << std::endl;
727: }
728: }
729: PetscSynchronizedPrintf(comm, txt.str().c_str());
730: PetscSynchronizedFlush(comm);
731: };
732: };
734: template<typename Bundle_, typename Sieve_>
735: class SupportSizeSection : public ALE::ParallelObject {
736: public:
737: typedef Bundle_ bundle_type;
738: typedef Sieve_ sieve_type;
739: typedef typename sieve_type::source_type point_type;
740: typedef typename sieve_type::target_type value_type;
741: protected:
742: Obj<bundle_type> _bundle;
743: Obj<sieve_type> _sieve;
744: value_type _size;
745: int _minDepth;
746: public:
747: 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) {};
748: virtual ~SupportSizeSection() {};
749: public:
750: bool hasPoint(const point_type& point) {return true;};
751: const value_type *restrictPoint(const point_type& p) {
752: if ((this->_minDepth == 0) || (this->_bundle->depth(p) >= this->_minDepth)) {
753: this->_size = this->_sieve->support(p)->size();
754: } else {
755: this->_size = 0;
756: }
757: return &this->_size;
758: };
759: public:
760: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
761: ostringstream txt;
762: int rank;
764: if (comm == MPI_COMM_NULL) {
765: comm = this->comm();
766: rank = this->commRank();
767: } else {
768: MPI_Comm_rank(comm, &rank);
769: }
770: if (name == "") {
771: if(rank == 0) {
772: txt << "viewing a SupportSizeSection" << std::endl;
773: }
774: } else {
775: if(rank == 0) {
776: txt << "viewing SupportSizeSection '" << name << "'" << std::endl;
777: }
778: }
779: PetscSynchronizedPrintf(comm, txt.str().c_str());
780: PetscSynchronizedFlush(comm);
781: };
782: };
784: template<typename Sieve_>
785: class SupportSection : public ALE::ParallelObject {
786: public:
787: typedef Sieve_ sieve_type;
788: typedef typename sieve_type::source_type point_type;
789: typedef typename sieve_type::target_type value_type;
790: typedef PartitionDomain<sieve_type> chart_type;
791: protected:
792: Obj<sieve_type> _sieve;
793: int _supportSize;
794: value_type *_support;
795: chart_type _domain;
796: void ensureSupport(const int size) {
797: if (size > this->_supportSize) {
798: if (this->_support) delete [] this->_support;
799: this->_supportSize = size;
800: this->_support = new value_type[this->_supportSize];
801: }
802: };
803: public:
804: SupportSection(const Obj<sieve_type>& sieve) : ParallelObject(MPI_COMM_SELF, sieve->debug()), _sieve(sieve), _supportSize(-1), _support(NULL) {};
805: virtual ~SupportSection() {if (this->_support) delete [] this->_support;};
806: public:
807: const chart_type& getChart() {return this->_domain;};
808: bool hasPoint(const point_type& point) {return true;};
809: const value_type *restrictPoint(const point_type& p) {
810: const Obj<typename sieve_type::traits::supportSequence>& support = this->_sieve->support(p);
811: int s = 0;
813: this->ensureSupport(support->size());
814: for(typename sieve_type::traits::supportSequence::iterator s_iter = support->begin(); s_iter != support->end(); ++s_iter) {
815: this->_support[s++] = *s_iter;
816: }
817: return this->_support;
818: };
819: public:
820: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
821: ostringstream txt;
822: int rank;
824: if (comm == MPI_COMM_NULL) {
825: comm = this->comm();
826: rank = this->commRank();
827: } else {
828: MPI_Comm_rank(comm, &rank);
829: }
830: if (name == "") {
831: if(rank == 0) {
832: txt << "viewing a SupportSection" << std::endl;
833: }
834: } else {
835: if(rank == 0) {
836: txt << "viewing SupportSection '" << name << "'" << std::endl;
837: }
838: }
839: PetscSynchronizedPrintf(comm, txt.str().c_str());
840: PetscSynchronizedFlush(comm);
841: };
842: };
844: template<typename Sieve_, typename Section_>
845: class ArrowSection : public ALE::ParallelObject {
846: public:
847: typedef Sieve_ sieve_type;
848: typedef Section_ section_type;
849: typedef typename sieve_type::target_type point_type;
850: typedef typename section_type::point_type arrow_type;
851: typedef typename section_type::value_type value_type;
852: protected:
853: Obj<sieve_type> _sieve;
854: Obj<section_type> _section;
855: int _coneSize;
856: value_type *_cone;
857: void ensureCone(const int size) {
858: if (size > this->_coneSize) {
859: if (this->_cone) delete [] this->_cone;
860: this->_coneSize = size;
861: this->_cone = new value_type[this->_coneSize];
862: }
863: };
864: public:
865: 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) {};
866: virtual ~ArrowSection() {if (this->_cone) delete [] this->_cone;};
867: public:
868: bool hasPoint(const point_type& point) {return this->_sieve->baseContains(point);};
869: const value_type *restrictPoint(const point_type& p) {
870: const Obj<typename sieve_type::traits::coneSequence>& cone = this->_sieve->cone(p);
871: int c = -1;
873: this->ensureCone(cone->size());
874: for(typename sieve_type::traits::coneSequence::iterator c_iter = cone->begin(); c_iter != cone->end(); ++c_iter) {
875: this->_cone[++c] = this->_section->restrictPoint(arrow_type(*c_iter, p))[0];
876: }
877: return this->_cone;
878: };
879: public:
880: void view(const std::string& name, MPI_Comm comm = MPI_COMM_NULL) const {
881: ostringstream txt;
882: int rank;
884: if (comm == MPI_COMM_NULL) {
885: comm = this->comm();
886: rank = this->commRank();
887: } else {
888: MPI_Comm_rank(comm, &rank);
889: }
890: if (name == "") {
891: if(rank == 0) {
892: txt << "viewing a ConeSection" << std::endl;
893: }
894: } else {
895: if(rank == 0) {
896: txt << "viewing ConeSection '" << name << "'" << std::endl;
897: }
898: }
899: PetscSynchronizedPrintf(comm, txt.str().c_str());
900: PetscSynchronizedFlush(comm);
901: };
902: };
903: }
904: }
905: #endif