Actual source code: SieveAlgorithms.hh
petsc-3.4.5 2014-06-29
1: #ifndef included_ALE_SieveAlgorithms_hh
2: #define included_ALE_SieveAlgorithms_hh
4: #ifndef included_ALE_Sieve_hh
5: #include <sieve/Sieve.hh>
6: #endif
8: namespace ALE {
9: template<typename Bundle_>
10: class SieveAlg {
11: public:
12: typedef Bundle_ bundle_type;
13: typedef typename bundle_type::sieve_type sieve_type;
14: typedef typename sieve_type::point_type point_type;
15: typedef typename sieve_type::coneSet coneSet;
16: typedef typename sieve_type::coneArray coneArray;
17: typedef typename sieve_type::supportSet supportSet;
18: typedef typename sieve_type::supportArray supportArray;
19: typedef typename bundle_type::arrow_section_type arrow_section_type;
20: typedef std::pair<point_type, int> oriented_point_type;
21: typedef ALE::array<oriented_point_type> orientedConeArray;
22: protected:
23: typedef MinimalArrow<point_type, point_type> arrow_type;
24: typedef std::pair<arrow_type, int> oriented_arrow_type;
25: typedef ALE::array<oriented_arrow_type> orientedArrowArray;
26: public:
27: static Obj<coneArray> closure(const Obj<bundle_type>& bundle, const point_type& p) {
28: return closure(bundle.ptr(), bundle->getArrowSection("orientation"), p);
29: };
30: static Obj<coneArray> closure(const Bundle_ *bundle, const Obj<arrow_section_type>& orientation, const point_type& p) {
31: const Obj<sieve_type>& sieve = bundle->getSieve();
32: const int depth = bundle->depth();
33: Obj<orientedArrowArray> cone = new orientedArrowArray();
34: Obj<orientedArrowArray> base = new orientedArrowArray();
35: Obj<coneArray> closure = new coneArray();
36: coneSet seen;
38: // Cone is guarateed to be ordered correctly
39: const Obj<typename sieve_type::traits::coneSequence>& initCone = sieve->cone(p);
41: closure->push_back(p);
42: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
43: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, p), 1));
44: closure->push_back(*c_iter);
45: }
46: for(int i = 1; i < depth; ++i) {
47: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
49: cone->clear();
50: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
51: const arrow_type& arrow = b_iter->first;
52: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source);
53: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0];
55: if (b_iter->second < 0) {
56: o = -(o+1);
57: }
58: if (o < 0) {
59: const int size = pCone->size();
61: if (o == -size) {
62: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
63: if (seen.find(*c_iter) == seen.end()) {
64: seen.insert(*c_iter);
65: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
66: closure->push_back(*c_iter);
67: }
68: }
69: } else {
70: const int numSkip = size + o;
71: int count = 0;
73: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
74: if (count < numSkip) continue;
75: if (seen.find(*c_iter) == seen.end()) {
76: seen.insert(*c_iter);
77: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
78: closure->push_back(*c_iter);
79: }
80: }
81: count = 0;
82: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
83: if (count >= numSkip) break;
84: if (seen.find(*c_iter) == seen.end()) {
85: seen.insert(*c_iter);
86: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
87: closure->push_back(*c_iter);
88: }
89: }
90: }
91: } else {
92: if (o == 1) {
93: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
94: if (seen.find(*c_iter) == seen.end()) {
95: seen.insert(*c_iter);
96: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
97: closure->push_back(*c_iter);
98: }
99: }
100: } else {
101: const int numSkip = o-1;
102: int count = 0;
104: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
105: if (count < numSkip) continue;
106: if (seen.find(*c_iter) == seen.end()) {
107: seen.insert(*c_iter);
108: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
109: closure->push_back(*c_iter);
110: }
111: }
112: count = 0;
113: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
114: if (count >= numSkip) break;
115: if (seen.find(*c_iter) == seen.end()) {
116: seen.insert(*c_iter);
117: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
118: closure->push_back(*c_iter);
119: }
120: }
121: }
122: }
123: }
124: }
125: return closure;
126: };
127: static Obj<orientedConeArray> orientedClosure(const Obj<bundle_type>& bundle, const point_type& p) {
128: return orientedClosure(bundle.ptr(), bundle->getArrowSection("orientation"), p);
129: };
130: static Obj<orientedConeArray> orientedClosure(const bundle_type *bundle, const Obj<arrow_section_type>& orientation, const point_type& p) {
131: const Obj<sieve_type>& sieve = bundle->getSieve();
132: const int depth = bundle->depth();
133: Obj<orientedArrowArray> cone = new orientedArrowArray();
134: Obj<orientedArrowArray> base = new orientedArrowArray();
135: Obj<orientedConeArray> closure = new orientedConeArray();
136: coneSet seen;
138: // Cone is guarateed to be ordered correctly
139: const Obj<typename sieve_type::traits::coneSequence>& initCone = sieve->cone(p);
141: closure->push_back(oriented_point_type(p, 0));
142: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
143: const arrow_type arrow(*c_iter, p);
145: cone->push_back(oriented_arrow_type(arrow, 1));
146: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(arrow)[0]));
147: }
148: for(int i = 1; i < depth; ++i) {
149: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
151: cone->clear();
152: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
153: const arrow_type& arrow = b_iter->first;
154: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source);
155: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0];
157: if (b_iter->second < 0) {
158: o = -(o+1);
159: }
160: if (o < 0) {
161: const int size = pCone->size();
163: if (o == -size) {
164: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
165: if (seen.find(*c_iter) == seen.end()) {
166: const arrow_type newArrow(*c_iter, arrow.source);
167: int pointO = orientation->restrictPoint(newArrow)[0];
169: seen.insert(*c_iter);
170: cone->push_back(oriented_arrow_type(newArrow, o));
171: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
172: }
173: }
174: } else {
175: const int numSkip = size + o;
176: int count = 0;
178: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
179: if (count < numSkip) continue;
180: if (seen.find(*c_iter) == seen.end()) {
181: const arrow_type newArrow(*c_iter, arrow.source);
182: int pointO = orientation->restrictPoint(newArrow)[0];
184: seen.insert(*c_iter);
185: cone->push_back(oriented_arrow_type(newArrow, o));
186: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
187: }
188: }
189: count = 0;
190: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
191: if (count >= numSkip) break;
192: if (seen.find(*c_iter) == seen.end()) {
193: const arrow_type newArrow(*c_iter, arrow.source);
194: int pointO = orientation->restrictPoint(newArrow)[0];
196: seen.insert(*c_iter);
197: cone->push_back(oriented_arrow_type(newArrow, o));
198: closure->push_back(oriented_point_type(*c_iter, pointO ? -(pointO+1): pointO));
199: }
200: }
201: }
202: } else {
203: if (o == 1) {
204: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
205: if (seen.find(*c_iter) == seen.end()) {
206: const arrow_type newArrow(*c_iter, arrow.source);
208: seen.insert(*c_iter);
209: cone->push_back(oriented_arrow_type(newArrow, o));
210: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
211: }
212: }
213: } else {
214: const int numSkip = o-1;
215: int count = 0;
217: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
218: if (count < numSkip) continue;
219: if (seen.find(*c_iter) == seen.end()) {
220: const arrow_type newArrow(*c_iter, arrow.source);
222: seen.insert(*c_iter);
223: cone->push_back(oriented_arrow_type(newArrow, o));
224: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
225: }
226: }
227: count = 0;
228: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
229: if (count >= numSkip) break;
230: if (seen.find(*c_iter) == seen.end()) {
231: const arrow_type newArrow(*c_iter, arrow.source);
233: seen.insert(*c_iter);
234: cone->push_back(oriented_arrow_type(newArrow, o));
235: closure->push_back(oriented_point_type(*c_iter, orientation->restrictPoint(newArrow)[0]));
236: }
237: }
238: }
239: }
240: }
241: }
242: return closure;
243: };
244: static Obj<coneArray> nCone(const Obj<bundle_type>& bundle, const point_type& p, const int n) {
245: const Obj<sieve_type>& sieve = bundle->getSieve();
246: const Obj<arrow_section_type>& orientation = bundle->getArrowSection("orientation");
247: const int height = std::min(n, bundle->height());
248: Obj<orientedArrowArray> cone = new orientedArrowArray();
249: Obj<orientedArrowArray> base = new orientedArrowArray();
250: Obj<coneArray> nCone = new coneArray();
251: coneSet seen;
253: if (height == 0) {
254: nCone->push_back(p);
255: return nCone;
256: }
258: // Cone is guarateed to be ordered correctly
259: const Obj<typename sieve_type::traits::coneSequence>& initCone = sieve->cone(p);
261: if (height == 1) {
262: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
263: nCone->push_back(*c_iter);
264: }
265: return nCone;
266: } else {
267: for(typename sieve_type::traits::coneSequence::iterator c_iter = initCone->begin(); c_iter != initCone->end(); ++c_iter) {
268: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, p), 1));
269: }
270: }
271: for(int i = 1; i < height; ++i) {
272: Obj<orientedArrowArray> tmp = cone; cone = base; base = tmp;
274: cone->clear();
275: for(typename orientedArrowArray::iterator b_iter = base->begin(); b_iter != base->end(); ++b_iter) {
276: const arrow_type& arrow = b_iter->first;
277: const Obj<typename sieve_type::traits::coneSequence>& pCone = sieve->cone(arrow.source);
278: typename arrow_section_type::value_type o = orientation->restrictPoint(arrow)[0];
280: if (b_iter->second < 0) {
281: o = -(o+1);
282: }
283: if (o < 0) {
284: const int size = pCone->size();
286: if (o == -size) {
287: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter) {
288: if (seen.find(*c_iter) == seen.end()) {
289: seen.insert(*c_iter);
290: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
291: if (i == height-1) nCone->push_back(*c_iter);
292: }
293: }
294: } else {
295: const int numSkip = size + o;
296: int count = 0;
298: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
299: if (count < numSkip) continue;
300: if (seen.find(*c_iter) == seen.end()) {
301: seen.insert(*c_iter);
302: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
303: if (i == height-1) nCone->push_back(*c_iter);
304: }
305: }
306: count = 0;
307: for(typename sieve_type::traits::coneSequence::reverse_iterator c_iter = pCone->rbegin(); c_iter != pCone->rend(); ++c_iter, ++count) {
308: if (count >= numSkip) break;
309: if (seen.find(*c_iter) == seen.end()) {
310: seen.insert(*c_iter);
311: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
312: if (i == height-1) nCone->push_back(*c_iter);
313: }
314: }
315: }
316: } else {
317: if (o == 1) {
318: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter) {
319: if (seen.find(*c_iter) == seen.end()) {
320: seen.insert(*c_iter);
321: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
322: if (i == height-1) nCone->push_back(*c_iter);
323: }
324: }
325: } else {
326: const int numSkip = o-1;
327: int count = 0;
329: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
330: if (count < numSkip) continue;
331: if (seen.find(*c_iter) == seen.end()) {
332: seen.insert(*c_iter);
333: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
334: if (i == height-1) nCone->push_back(*c_iter);
335: }
336: }
337: count = 0;
338: for(typename sieve_type::traits::coneSequence::iterator c_iter = pCone->begin(); c_iter != pCone->end(); ++c_iter, ++count) {
339: if (count >= numSkip) break;
340: if (seen.find(*c_iter) == seen.end()) {
341: seen.insert(*c_iter);
342: cone->push_back(oriented_arrow_type(arrow_type(*c_iter, arrow.source), o));
343: if (i == height-1) nCone->push_back(*c_iter);
344: }
345: }
346: }
347: }
348: }
349: }
350: return nCone;
351: };
352: static Obj<supportArray> star(const Obj<bundle_type>& bundle, const point_type& p) {
353: typedef MinimalArrow<point_type, point_type> arrow_type;
354: typedef typename ALE::array<arrow_type> arrowArray;
355: const Obj<sieve_type>& sieve = bundle->getSieve();
356: const Obj<arrow_section_type>& orientation = bundle->getArrowSection("orientation");
357: const int height = bundle->height();
358: Obj<arrowArray> support = new arrowArray();
359: Obj<arrowArray> cap = new arrowArray();
360: Obj<supportArray> star = new supportArray();
361: supportSet seen;
363: // Support is guarateed to be ordered correctly
364: const Obj<typename sieve_type::traits::supportSequence>& initSupport = sieve->support(p);
366: star->push_back(p);
367: for(typename sieve_type::traits::supportSequence::iterator s_iter = initSupport->begin(); s_iter != initSupport->end(); ++s_iter) {
368: support->push_back(arrow_type(p, *s_iter));
369: star->push_back(*s_iter);
370: }
371: for(int i = 1; i < height; ++i) {
372: Obj<arrowArray> tmp = support; support = cap; cap = tmp;
374: support->clear();
375: for(typename arrowArray::iterator b_iter = cap->begin(); b_iter != cap->end(); ++b_iter) {
376: const Obj<typename sieve_type::traits::supportSequence>& pSupport = sieve->support(b_iter->target);
377: const typename arrow_section_type::value_type o = orientation->restrictPoint(*b_iter)[0];
379: if (o == -1) {
380: for(typename sieve_type::traits::supportSequence::reverse_iterator s_iter = pSupport->rbegin(); s_iter != pSupport->rend(); ++s_iter) {
381: if (seen.find(*s_iter) == seen.end()) {
382: seen.insert(*s_iter);
383: support->push_back(arrow_type(b_iter->target, *s_iter));
384: star->push_back(*s_iter);
385: }
386: }
387: } else {
388: for(typename sieve_type::traits::supportSequence::iterator s_iter = pSupport->begin(); s_iter != pSupport->end(); ++s_iter) {
389: if (seen.find(*s_iter) == seen.end()) {
390: seen.insert(*s_iter);
391: support->push_back(arrow_type(b_iter->target, *s_iter));
392: star->push_back(*s_iter);
393: }
394: }
395: }
396: }
397: }
398: return star;
399: };
400: static Obj<supportArray> nSupport(const Obj<bundle_type>& bundle, const point_type& p, const int n) {
401: typedef MinimalArrow<point_type, point_type> arrow_type;
402: typedef typename ALE::array<arrow_type> arrowArray;
403: const Obj<sieve_type>& sieve = bundle->getSieve();
404: const Obj<arrow_section_type>& orientation = bundle->getArrowSection("orientation");
405: const int depth = std::min(n, bundle->depth());
406: Obj<arrowArray> support = new arrowArray();
407: Obj<arrowArray> cap = new arrowArray();
408: Obj<coneArray> nSupport = new supportArray();
409: supportSet seen;
411: if (depth == 0) {
412: nSupport->push_back(p);
413: return nSupport;
414: }
416: // Cone is guarateed to be ordered correctly
417: const Obj<typename sieve_type::traits::supportSequence>& initSupport = sieve->support(p);
419: if (depth == 1) {
420: for(typename sieve_type::traits::supportSequence::iterator s_iter = initSupport->begin(); s_iter != initSupport->end(); ++s_iter) {
421: nSupport->push_back(*s_iter);
422: }
423: return nSupport;
424: } else {
425: for(typename sieve_type::traits::supportSequence::iterator s_iter = initSupport->begin(); s_iter != initSupport->end(); ++s_iter) {
426: support->push_back(arrow_type(*s_iter, p));
427: }
428: }
429: for(int i = 1; i < depth; ++i) {
430: Obj<arrowArray> tmp = support; support = cap; cap = tmp;
432: support->clear();
433: for(typename arrowArray::iterator c_iter = cap->begin(); c_iter != cap->end(); ++c_iter) {
434: const Obj<typename sieve_type::traits::supportSequence>& pSupport = sieve->support(c_iter->source);
435: const typename arrow_section_type::value_type o = orientation->restrictPoint(*c_iter)[0];
437: if (o == -1) {
438: for(typename sieve_type::traits::supportSequence::reverse_iterator s_iter = pSupport->rbegin(); s_iter != pSupport->rend(); ++s_iter) {
439: if (seen.find(*s_iter) == seen.end()) {
440: seen.insert(*s_iter);
441: support->push_back(arrow_type(*s_iter, c_iter->source));
442: if (i == depth-1) nSupport->push_back(*s_iter);
443: }
444: }
445: } else {
446: for(typename sieve_type::traits::supportSequence::iterator s_iter = pSupport->begin(); s_iter != pSupport->end(); ++s_iter) {
447: if (seen.find(*s_iter) == seen.end()) {
448: seen.insert(*s_iter);
449: support->push_back(arrow_type(*s_iter, c_iter->source));
450: if (i == depth-1) nSupport->push_back(*s_iter);
451: }
452: }
453: }
454: }
455: }
456: return nSupport;
457: };
458: static Obj<sieve_type> Link(const Obj<bundle_type>& bundle, const point_type& p) {
459: const Obj<sieve_type>& link_sieve = new sieve_type(bundle->comm(), bundle->debug());
460: const Obj<sieve_type>& sieve = bundle->getSieve();
461: const int depth = bundle->depth(p);
462: const int interpolation_depth = bundle->height(p)+depth;
463: Obj<coneArray> nSupport = new supportArray();
464: supportSet seen;
465: //in either case, copy the closure of the cells surrounding the point to the new sieve
466: static Obj<supportArray> neighboring_cells = sieve->nSupport(sieve->nCone(p, depth), interpolation_depth);
467: static typename supportArray::iterator nc_iter = neighboring_cells->begin();
468: static typename supportArray::iterator nc_iter_end = neighboring_cells->end();
469: while (nc_iter != nc_iter_end) {
470: addClosure(sieve, link_sieve, *nc_iter);
471: nc_iter++;
472: }
473: if (interpolation_depth == 1) { //noninterpolated case
474: //remove the point, allowing the copied closure to contract to a surface.
476: if (depth != 0) {
477: static Obj<coneArray> point_cone = sieve->cone(p);
478: static typename coneArray::iterator pc_iter = point_cone->begin();
479: static typename coneArray::iterator pc_iter_end = point_cone->end();
480: while (pc_iter != pc_iter_end) {
481: link_sieve->removePoint(*pc_iter);
482: pc_iter++;
483: }
484: }
485: link_sieve->removePoint(p);
487: } else { //interpolated case: remove the point, its closure, and the support of that closure,
489: static Obj<supportArray> surrounding_support = sieve->Star(sieve->nCone(p, depth));
490: static typename supportArray::iterator ss_iter = surrounding_support->begin();
491: static typename supportArray::iterator ss_iter_end = surrounding_support->end();
492: while (ss_iter != ss_iter_end) {
493: link_sieve->removePoint(*ss_iter);
494: ss_iter++;
495: }
496: link_sieve->removePoint(p);
497: }
498: link_sieve->stratify();
499: return link_sieve;
500: };
501: };
502: }
504: #endif