Actual source code: SieveAlgorithms.hh
petsc-3.4.2 2013-07-02

      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