source: issm/trunk-jpl/externalpackages/petsc-dev/src/include/sieve/Sections.hh@ 11896

Last change on this file since 11896 was 11896, checked in by habbalf, 13 years ago

petsc-dev : Petsc development code in external packages. Mercurial updates
File size: 38.2 KB
Line 
1#ifndef included_ALE_Sections_hh
2#define included_ALE_Sections_hh
3
4namespace 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 };
39
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 };
78
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;
105
106 if (size > maxSize) {
107 const value_type dummy(0);
108
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;
135
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 };
142
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 };
179
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 };
218
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 };
264
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 };
318
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 };
355
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 };
397
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 };
436
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 };
464
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;
482
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);
490
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);
494
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;
516
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;
535
536 txt << "[" << this->commRank() << "]: Partition " << partition << " size " << size << std::endl;
537 }
538 PetscSynchronizedPrintf(comm, txt.str().c_str());
539 PetscSynchronizedFlush(comm);
540 };
541 };
542
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 };
553
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;
575
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);
583
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);
587
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;
625
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;
644
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 };
654
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;
695
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 };
715
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;
744
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;
755
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 };
775
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;
805
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 };
825
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;
854
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;
865
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 };
885
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;
914
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;
925
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
Note: See TracBrowser for help on using the repository browser.